[CP2K:4532] Re: Comparing the energy from cp2k and gromacs

[Fahimeh Baftizadeh]

Dear Peter,

Thanks for the useful reply. I really appreciate it.

In order to be sure, I took a sample crd and top file from cp2k sample 
files. They are called ace_ala_nme.crd and ace_ala_nme.top  
All the charges and force field parameters are correct while loaded in 
cp2k. I checked them carefully. 

Then I preformed a single step energy calculation with cp2k, the run_type 
is ENERGY ... so it doesn't move anything. just compute the energy. 

for gromacs, I used the same top and crd file, just convert them with 
amb2gmx tool to gromacs format. This is very common thing in the community. 
Then I performed a rerun on the structure. and I compute the energy of that 
crd file. 

As you correctly mentioned I thought that  Angle and Torsion energies 
should match perfectly between gromacs and cp2k. But they dont match. 

These are the values in cp2k :  ANGLE  = 0.3620 (kcal/mol)   TORSION = 
 8.1071 (kcal/mol)
These are the values from gromacs:  ANGLE= 0.442935141 (kcal/mol) 
TORSION=8.11232826 (kcal/mol)

As you see bonded terms don't match perfectly. now lets go to the other 
terms, in both programs i am using pme:
In gromacs there are short range and long range lennard jones terms and 
long range and short range coulomb terms, while in cp2k the notation is 
different. however the total energy of the molecule in two packages differs 
as the following:

ENERGY cp2k = -0.102733788208951 Hartree = -269.72 (kj/mol)
ENERGY gromacs= -158.67 (kj/mol)

What do you think? this doesn't seem a good sign! I am not doing anything 
complicated! it is very trivial!
Am i missing something?

Fahimeh


On Wednesday, August 7, 2013 6:02:13 AM UTC-4, Peter wrote:
>
> Hello Fahimeh,
>
> Check the energy components: bonds, angles, torsion angles, coulomb and 
> vdw. I believe the bonded interactions energies should agree perfectly, 
> otherwise there is a mismatch between topologies and/or forcefield you use 
> in gromacs and cp2k. If those are ok, then I guess this discrepancy is due 
> to differences in PME electrostatics calculation/implementation.
>
> Regards,
> Peter
>
>
> On Tue, Aug 6, 2013 at 8:00 PM, Fahimeh Baftizadeh <
> fahime... at googlemail.com <javascript:>> wrote:
>
>> Hello everyone, 
>>
>> I thought maybe it will be helpful if I just put here my input files and 
>> someone else run it if it is possible for you. In fact it seems very 
>> trivial but I checked the parameters and I don't understand why the Energy 
>> is positive! 
>>
>> I am doing just a single step energy calculation of a molecule. 
>> Thanks for your help
>>
>> Fahimeh
>>
>>
>>
>> On Monday, August 5, 2013 4:53:08 PM UTC-4, Fahimeh Baftizadeh wrote:
>>>
>>> Hello,
>>>
>>> I have a system which contains only one molecule in gas phase. I wanna 
>>> compute its energy with cp2k and gromacs (for a test).
>>> I am reading a pdb file and the force filed is in amber format, this is 
>>> cp2k input file:
>>>
>>> &FORCE_EVAL
>>>   METHOD FIST
>>>   &MM
>>>     &FORCEFIELD
>>>       &SPLINE
>>>         EMAX_ACCURACY 500.0
>>>         EMAX_SPLINE    1000.0
>>>         EPS_SPLINE [hartree] 1.00000000E-07
>>>       &END
>>>       parm_file_name /home/fbafti/work/nucleation/**
>>> test-ff/single-mol/amber-1res.**top
>>>       parmtype AMBER
>>>     &END FORCEFIELD
>>>     &POISSON
>>>      &EWALD
>>>        EWALD_TYPE pme
>>>        RCUT [angstrom] 8.0
>>>        EPSILON 1E-6
>>>      &END EWALD
>>>     &END POISSON
>>>     &PRINT
>>>       &FF_INFO
>>>       &END
>>>     &END
>>>   &END MM
>>>   &SUBSYS
>>>     &CELL
>>>       ABC 125.990  122.910  129.640
>>>     &END
>>>     &TOPOLOGY
>>>      COORD_FILE_NAME /home/fbafti/work/nucleation/**
>>> test-ff/single-mol/non-per.pdb
>>>      COORDINATE PDB
>>>      CONN_FILE_NAME  /home/fbafti/work/nucleation/**
>>> test-ff/single-mol/amber-1res.**top
>>>      CONN_FILE_FORMAT  AMBER
>>>     &END TOPOLOGY
>>>   &END SUBSYS
>>>   STRESS_TENSOR ANALYTICAL
>>> &END FORCE_EVAL
>>> &GLOBAL
>>>   PROJECT paracetamol_single_mol
>>>   RUN_TYPE ENERGY
>>> &END GLOBAL
>>>
>>>
>>> I put a very large box size, just to be sure that there are no periodic 
>>> images in the energy estimation. The energy output of CP2K has a value 
>>> which is larger than zero and it is strange.
>>> I performed a single step steepest decent with gromacs and I got another 
>>> value which is minus and makes more sense.
>>> I already checked the electrostatic parameters and vdw ... I think they 
>>> are all similar while the energy difference coming from gromacs and 
>>> cp2k differs by almost 500 kj/mol and in gromacs is minus and in cp2k is a 
>>> plus value.
>>>
>>> This is the gromacs input file: 
>>>
>>> integrator               = steep
>>> comm_mode                = Linear
>>> dt                       = 0.001
>>> nsteps                   = 1
>>> ; Output frequency for coords (x), velocities (v) and forces (f) =
>>> nstxout                  = 1
>>> nstvout                  = 1
>>> ; Output frequency for energies to log file and energy file =
>>> nstlog                   = 1
>>> nstenergy                = 1
>>> nstxtcout                = 1
>>> xtc-precision            = 1000000
>>> xtc_grps                 = System
>>> energygrps               = System
>>> pbc                      = xyz
>>> nstlist                  = 10
>>> epsilon_r                = 1.
>>> ns_type                  = grid ; or grid I don't know
>>> coulombtype              = pme
>>> vdwtype                  = Cut-off
>>> fourierspacing           = 0.12
>>> ; EWALD/PME/PPPM parameters
>>> pme_order                = 4
>>> ewald_rtol               = 1e-05
>>> epsilon_surface          = 0
>>> optimize_fft             = yes
>>> rlist                    = 0.8
>>> rcoulomb                 = 0.8
>>> rvdw                     = 0.8
>>> emtol                    = 1.0
>>>
>>> Can anybody help me with that?
>>> Thanks
>>> Fahimeh
>>>
>>  -- 
>> You received this message because you are subscribed to the Google Groups 
>> "cp2k" group.
>> To unsubscribe from this group and stop receiving emails from it, send an 
>> email to cp2k+... at googlegroups.com <javascript:>.
>> To post to this group, send email to cp... at googlegroups.com <javascript:>
>> .
>> Visit this group at http://groups.google.com/group/cp2k.
>> For more options, visit https://groups.google.com/groups/opt_out.
>>  
>>  
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20130807/da2b073b/attachment.htm>