Skip to content

[Codescan] Release or reuse EDA-FF parameter arrays before re-entering the analysis #34

Description

@Stardust0831

Global scan provenance

This issue was found by the Codex global repository scan of every non-unit-test file in Stardust0831/Multiwfn, pinned to commit d2b770a531eb7039153198d04520b5992f15c49d.

Problem

parmA and parmB are allocatable variables in the saved EDA_FF_mod module state:

Multiwfn/EDA.f90

Lines 44 to 48 in d2b770a

!------ Module for EDA-FF
module EDA_FF_mod
integer :: ielemode=1,ivdwmode=2
real*8,allocatable :: parmA(:),parmB(:)
end module

Every entry into EDA_forcefield() allocates both arrays unconditionally:

Multiwfn/EDA.f90

Lines 50 to 66 in d2b770a

subroutine EDA_forcefield
use defvar
use util
use EDA_FF_mod
implicit real*8 (a-h,o-z)
integer :: nfrag=0,ishowatmpair=0,ioutatmpqr=0
character c200tmp*200,c200tmp2*200,c2000tmp*2000
integer,allocatable :: frag(:,:),fragnatm(:)
real*8,allocatable :: elemat(:,:),repmat(:,:),dispmat(:,:),totmat(:,:) !Interfragment interaction matrix of electrostatic, repulsive, dispersion and total interaction
real*8 eleatm(ncenter),repatm(ncenter),dispatm(ncenter),totatm(ncenter)
!vdW parameter of each atom, for UFF they are well depths and vdW distance; For AMBER/GAFF they are well depths and vdW radii
character(len=2) c2tmp,FFtype(ncenter) !Force field atom types, only needed by AMBER and GAFF
FFtype="?"
if (ifiletype/=4) a%charge=0 !If the input file is not chg or pqr format, assume atomic charges to be zero
allocate(parmA(ncenter),parmB(ncenter))

Selecting 0 Return exits the subroutine without deallocating them:

Multiwfn/EDA.f90

Lines 81 to 122 in d2b770a

write(*,*) "0 Return"
write(*,*) "1 Start analysis"
if (allocated(frag)) then
write(*,"(a,i4,a)") " 2 Redefine fragments, current:",nfrag," fragments"
else
write(*,"(a)") " 2 Define fragments (undefined currently)"
end if
if (ivdwmode==1) then !UFF
write(*,*) "3 Load atomic charges for current system"
else !AMBER&GAFF
write(*,*) "3 Load atomic types and charges for current system"
end if
if (ivdwmode==1) then !UFF
write(*,*) "4 Show current atomic charges"
else !AMBER&GAFF
write(*,*) "4 Show current atomic types and charges"
end if
read(*,*) isel
if (isel==-4) then
if (ioutatmpqr==1) then
ioutatmpqr=0
else
ioutatmpqr=1
end if
else if (isel==-3) then
if (ishowatmpair==1) then
ishowatmpair=0
else
ishowatmpair=1
end if
else if (isel==-2) then
write(*,*) "1: 1/r interaction potential"
write(*,*) "2: 1/r^2 interaction potential"
read(*,*) ielemode
else if (isel==-1) then
write(*,*) "1: UFF van der Waals model"
write(*,*) "2: AMBER99 & GAFF van der Waals model"
read(*,*) ivdwmode
else if (isel==0) then
return

The parent EDA menu remains active and allows option 1 to call EDA_forcefield() again in the same process:

Multiwfn/EDA.f90

Lines 1 to 30 in d2b770a

!-------- Main interface of various energy decomposition analyses
subroutine EDA_main
use defvar
implicit real*8 (a-h,o-z)
do while(.true.)
write(*,*)
write(*,*) " ============ Energy decomposition analysis ============ "
write(*,*) "0 Return"
write(*,*) "1 Energy decomposition analysis based on molecular forcefield (EDA-FF)"
write(*,*) "2 Shubin Liu's energy decomposition analysis (Gaussian is needed)"
write(*,*) "3 sobEDA and sobEDAw energy decomposition analyses (J. Phys. Chem. A, 127, 7023 (2023))"
write(*,*) "4 Analysis of atomic contribution to dispersion energy"
! write(*,*) "2 Mayer energy decomposition analysis"
! write(*,*) "3 Fuzzy space based energy decomposition analysis"
read(*,*) infuncsel2
if (infuncsel2==0) then
return
else if (infuncsel2==1) then
call EDA_forcefield
else if (infuncsel2==2) then
call EDA_SBL
else if (infuncsel2==3) then
write(*,"(a)") " This kind of analysis needs using shell script. Please check detailed sobEDA/sobEDAw tutorial: http://sobereva.com/soft/sobEDA_tutorial.zip"
write(*,*) "Also see original paper: J. Phys. Chem. A, 127, 7023 (2023) DOI: 10.1021/acs.jpca.3c04374"
write(*,*) "Press ENTER to return"
read(*,*)
else if (infuncsel2==4) then
call atomdispcontri
end if
end do

The second entry therefore executes allocate(parmA(...), parmB(...)) while both allocatables are already allocated, which is a Fortran runtime error.

Impact

A normal interactive sequence, EDA -> EDA-FF -> Return -> EDA-FF, can terminate the program instead of reopening the analysis. Reloading a structure with a different number of atoms also leaves arrays with stale dimensions until the failing allocation is addressed.

Suggested direction

  • Prefer local allocatables if these parameters do not need to survive the subroutine.
  • Otherwise, deallocate them on every return path or guard/reallocate when their size differs from ncenter.
  • Add a cleanup routine for EDA-FF module state and call it from whole-system cleanup.
  • Add a functional test that enters and leaves EDA-FF twice in one process.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions