[Freeon-users] Bug in QCTC.f90

[Jose R. Valverde]

Dear FreeON developers,

        I keep trying to use the core Hamiltonian guess (it's mainly a
"sports" issue, I want to be able to compare results with other codes)
and, yes, I know it is nit-picking, sorry for being such a nuisance.

        I just discovered a minor bug in QCTC.f90. The problem appears
when trying to use the core guess and results in the error 
"Counting error in IntegrateNCollate".

        I added a couple of "write(*,*)" and found that the value
Collate (Density.f90) gets is absurd and varies from run to run. 
Looking into QCTC.f90, just before the call to Collate, NLink is
added a value from GM%NAtms, but since for the core guess electrons
are not used, NLink seems (to me) to be used uninitialized here and
at least with gfortran, it may contain any (invalid) value.

        BTW, in Density.f90 one of the errors is "To many" instead
of "Too many.." but that's harmless.

        Adding a statement "NLink=0" at the beginning of QCTC.f90
just after the variable declarations fixes this issue for QCTC.
Using the HeH+-force.inp file:

        The calculation reaches to the energies and final KinE, Etot,
Exch energies obtained match TOTAL KINETIC ENERGY, TOTAL ENERGY and
ELECTRON-ELECTRON POTENTIAL ENERGY/TWO ELECTRON ENERGY and the sum
of E_el_tot, E_Nuc_tot and E_es_tot gives me approx the
NUCLEUS-ELECTRON POTENTIAL ENERGY (-1.0031005854504621D+01 vs.
-10.2932685839).

        That's about all the coincidences I can get, but at least
the main energies are OK.


        Now, for JForce, the problem is different, as NLink is
correctly initialized to 0 in MakeRhhoList which is called before
NLink is added GM%NAtms in JForce.f90. Comparing what Collate
expects, which is the exact value of GM%NAtms, and the value of
NLink before it is added GM%NAtms, I suppose that with the
core guess, either MakeRhoList should not be called, or either
GM%NAtms should be assigned and not added to NLink.

        If I change in JForce.f90 and use HeH+-force.inp as reference

  CALL MakeRhoList(GM,BS,P,NLink,RhoHead,'JForce',NoWrap_O=NoWrap)
  ! Add in the nuclear charges
  CALL AddNukes(GM,RhoHead)
  NLink=NLink+GM%NAtms
  ! Load density into arrays and delete the linked list
  CALL Collate(GM,RhoHead,Rho,'JForce',RhoPoles,NLink)

by

  CALL MakeRhoList(GM,BS,P,NLink,RhoHead,'JForce',NoWrap_O=NoWrap)
  ! Add in the nuclear charges
  CALL AddNukes(GM,RhoHead)
+  IF(SCFActn=='GuessEqCore')THEN
+    NLink=GM%NAtms
+  ELSE
    NLink=NLink+GM%NAtms
+  ENDIF
  ! Load density into arrays and delete the linked list
  CALL Collate(GM,RhoHead,Rho,'JForce',RhoPoles,NLink)

or alternately to

+  IF(SCFActn=='GuessEqCore')THEN
     CALL MakeRhoList(GM,BS,P,NLink,RhoHead,'JForce',NoWrap_O=NoWrap)
+  ELSE
+     NLink=0
+  ENDIF
  ! Add in the nuclear charges
  CALL AddNukes(GM,RhoHead)
  NLink=NLink+GM%NAtms
  ! Load density into arrays and delete the linked list
  CALL Collate(GM,RhoHead,Rho,'JForce',RhoPoles,NLink)


        Then the calculation completes, but the forces reported do
not seem to match any gradient value reported by GAMESS at all. In
both cases I seem to get the same results. So, I am a bit lost here,
I don't know if these values are OK and how to relate them to GAMESS
output:

From FreeON:
Force                            :: Force
Force                            :: Atom  Z  Forces (eV/A)
Force                            :: 1 2 -.9809040287581499D+02 
-.0000000000000000D+00 -.0000000000000000D+00
Force                            :: 2 1 0.1208438202044161D+03 
-.0000000000000000D+00 -.0000000000000000D+00
Force : Clone 1                  :: Positions (in A) and forces (in eV/A), |F| 
= 0.15564368D+03 eV/A, RMSd(F) = 0.11005670D+03 eV/A, max(F_{i}) = 
0.12084382D+03 eV/A
Force : Clone 1                  ::  HE    0.0000000000000000     
0.0000000000000000     0.0000000000000000    -98.0904028758149877    
-0.0000000000000000    -0.0000000000000000
Force : Clone 1                  ::  H     0.3704240460129999   
0.0000000000000000     0.0000000000000000  120.8438202044161045   
-0.0000000000000000    -0.0000000000000000

From GAMESS:
                         ----------------------
                         GRADIENT OF THE ENERGY
                         ----------------------

      ATOM                 E'X               E'Y               E'Z 
    1 HELIUM           2.001565573       0.000000000       0.000000000
    2 HYDROGEN        -2.001565573       0.000000000       0.000000000

                   MAXIMUM GRADIENT =    2.001565573
                       RMS GRADIENT =    1.155604422
  $VIB   
          IVIB=   0 IATOM=   0 ICOORD=   0 E=       -2.5168354952
  2.001565573E+00 0.000000000E+00 0.000000000E+00-2.001565573E+00 
0.000000000E+00
  0.000000000E+00
  4.434554816E-01 0.000000000E+00 0.000000000E+00


        Any suggestions? Where am I wrong here?

                                j
-- 
                        EMBnet/CNB
                Scientific Computing Service
        Solving all your computer needs for Scientific
                        Research.

                http://bioportal.cnb.csic.es
                  http://www.es.embnet.org