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

[Peter Mamonov]

Discrepancy in angles energies is quite strange: when I converted gromacs
TPRs to PSF format and fed them to cp2k, bonds and angles energies matched
perfectly (however there were minor differences in torsions energy). To
deal with this problem you should probably try to write those topologies/ff
by hand for some small system and compare the results, since I guess there
might be a problem with the conversion script.

Concerning PME energies, I guess you won't achieve a match between absolute
values due to differences in implementations. However what you really
should be concerned about are energy gradients (forces) produced by your
models. So try to compare forces calculated by gromacs and cp2k and if they
match to a reasonable degree then you are done )

Peter


On Wed, Aug 7, 2013 at 8:38 PM, Fahimeh Baftizadeh <
fahimeh.b... at googlemail.com> wrote:

> 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 <fahimeh.b...@**
>> googlemail.com> 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/**t**
>>>> est-ff/single-mol/amber-1res.**t**op
>>>>       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/**t**
>>>> est-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+uns... at googlegroups.**com.
>>> To post to this group, send email to cp... at googlegroups.com.
>>>
>>> Visit this group at http://groups.google.com/**group/cp2k<http://groups.google.com/group/cp2k>
>>> .
>>> For more options, visit https://groups.google.com/**groups/opt_out<https://groups.google.com/groups/opt_out>
>>> .
>>>
>>>
>>>
>>
>>  --
> 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+uns... at googlegroups.com.
> To post to this group, send email to cp... at googlegroups.com.
> 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/063a3dab/attachment.htm>