[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Freeon-users] Bug in QCTC.f90
From: |
Jose R. Valverde |
Subject: |
[Freeon-users] Bug in QCTC.f90 |
Date: |
Thu, 1 Mar 2012 17:01:14 +0100 |
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Freeon-users] Bug in QCTC.f90,
Jose R. Valverde <=