Comparing the energy from cp2k and gromacs

[Fahimeh Baftizadeh]

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

[Fahimeh Baftizadeh]

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

[Peter Mamonov]

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

--
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...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at http://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[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

[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