<div dir="ltr">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.<div>
<br></div><div>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 )</div>
<div><br></div><div>Peter<br><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Aug 7, 2013 at 8:38 PM, Fahimeh Baftizadeh <span dir="ltr"><<a href="mailto:fahimeh.b...@googlemail.com" target="_blank">fahimeh.b...@googlemail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Dear Peter,<div><br></div><div>Thanks for the useful reply. I really appreciate it.</div>
<div><br></div><div>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 <span style="font-size:13px"> </span></div><div>All the charges and force field parameters are correct while loaded in cp2k. I checked them carefully. </div>
<div><br></div><div>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. </div><div><br></div><div>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. </div>
<div><br></div><div>As you correctly mentioned I thought that <span style="background-color:rgb(255,153,0)"> Angle and Torsion</span> energies should match perfectly between gromacs and cp2k. But they dont match. </div><div>
<br></div><div>These are the values in cp2k : ANGLE = 0.3620 (kcal/mol) <span style="font-size:13px">TORSION = 8.1071 (kcal/mol)</span></div><div>These are the values from gromacs: ANGLE= 0.442935141 (kcal/mol) TORSION=8.11232826 (kcal/mol)</div>
<div><br></div><div>As you see bonded terms don't match perfectly. now lets go to the other terms, in both programs i am using pme:</div><div>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:</div>
<div><br></div><div>ENERGY cp2k = -0.102733788208951 Hartree = -269.72 (kj/mol)</div><div>ENERGY gromacs= -158.67 (kj/mol)</div><div><br></div><div>What do you think? this doesn't seem a good sign! I am not doing anything complicated! it is very trivial!</div>
<div>Am i missing something?</div><div><br></div><div>Fahimeh</div><div><br></div><div><div class="im"><br>On Wednesday, August 7, 2013 6:02:13 AM UTC-4, Peter wrote:</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div class="im"><div dir="ltr"><div>Hello Fahimeh,</div><div><br></div>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.<div>
<br></div><div>Regards,</div><div>Peter</div></div></div><div><br><br><div class="gmail_quote"><div><div class="h5">On Tue, Aug 6, 2013 at 8:00 PM, Fahimeh Baftizadeh <span dir="ltr"><<a>fahimeh.b...@<u></u>googlemail.com</a>></span> wrote:<br>
</div></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div><div class="h5">Hello everyone, <div><br>
</div><div>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! </div>
<div><br></div><div>I am doing just a single step energy calculation of a molecule. </div><div>Thanks for your help</div><span><font color="#888888"><div><br></div><div>Fahimeh</div></font></span></div></div><div>
<div><div><div class="h5"><div><br></div><div><br><br>On Monday, August 5, 2013 4:53:08 PM UTC-4, Fahimeh Baftizadeh wrote:<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Hello,<div><br></div><div>I have a system which contains only one molecule in gas phase. I wanna compute its energy with cp2k and gromacs (for a test).</div><div>I am reading a pdb file and the force filed is in amber format, this is cp2k input file:</div>
<div><br></div><div><div>&FORCE_EVAL</div><div> METHOD FIST</div><div> &MM</div><div> &FORCEFIELD</div><div> &SPLINE</div><div> EMAX_ACCURACY 500.0</div><div> EMAX_SPLINE 1000.0</div>
<div> EPS_SPLINE [hartree] 1.00000000E-07</div><div> &END</div><div> parm_file_name /home/fbafti/work/nucleation/<u></u>t<u></u>est-ff/single-mol/amber-1res.<u></u>t<u></u>op</div><div> parmtype AMBER</div>
<div>
&END FORCEFIELD</div><div> &POISSON</div><div> &EWALD</div><div> EWALD_TYPE pme</div><div> RCUT [angstrom] 8.0</div><div> EPSILON 1E-6</div><div> &END EWALD</div><div> &END POISSON</div>
<div> &PRINT</div><div> &FF_INFO</div><div> &END</div><div> &END</div><div> &END MM</div><div> &SUBSYS</div><div> &CELL</div><div> ABC 125.990 122.910 129.640</div>
<div>
&END</div><div> &TOPOLOGY</div><div> COORD_FILE_NAME /home/fbafti/work/nucleation/<u></u>t<u></u>est-ff/single-mol/non-per.pdb</div><div> COORDINATE PDB</div><div> CONN_FILE_NAME /home/fbafti/work/nucleation/<u></u><u></u>test-ff/single-mol/amber-1res.<u></u><u></u>top</div>
<div> CONN_FILE_FORMAT AMBER</div><div> &END TOPOLOGY</div><div> &END SUBSYS</div><div> STRESS_TENSOR ANALYTICAL</div><div>&END FORCE_EVAL</div><div>&GLOBAL</div><div> PROJECT paracetamol_single_mol</div>
<div> RUN_TYPE ENERGY</div><div>&END GLOBAL</div></div><div><br></div><div><br></div><div>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.</div>
<div>I performed a single step steepest decent with gromacs and I got another value which is minus and makes more sense.</div><div>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.</div>
<div><br></div><div>This is the gromacs input file: </div><div><div><br></div><div>integrator = steep</div><div>comm_mode = Linear</div><div>dt = 0.001</div><div>nsteps = 1</div>
<div>; Output frequency for coords (x), velocities (v) and forces (f) =</div><div>nstxout = 1</div><div>nstvout = 1</div><div>; Output frequency for energies to log file and energy file =</div>
<div>nstlog = 1</div><div>nstenergy = 1</div><div>nstxtcout = 1</div><div>xtc-precision = 1000000</div><div>xtc_grps = System</div><div>energygrps = System</div>
<div>pbc = xyz</div><div>nstlist = 10</div><div>epsilon_r = 1.</div><div>ns_type = grid ; or grid I don't know</div><div>coulombtype = pme</div>
<div>vdwtype = Cut-off</div><div>fourierspacing = 0.12</div><div>; EWALD/PME/PPPM parameters</div><div>pme_order = 4</div><div>ewald_rtol = 1e-05</div><div>epsilon_surface = 0</div>
<div>optimize_fft = yes</div><div>rlist = 0.8</div><div>rcoulomb = 0.8</div><div>rvdw = 0.8</div><div>emtol = 1.0</div></div><div><br>
</div>
<div>Can anybody help me with that?</div><div>Thanks</div><div>Fahimeh</div></blockquote></div>
<p></p>
-- <br>
You received this message because you are subscribed to the Google Groups "cp2k" group.<br></div></div>
To unsubscribe from this group and stop receiving emails from it, send an email to <a>cp2k+uns...@googlegroups.<u></u>com</a>.<br>
To post to this group, send email to <a>cp...@googlegroups.com</a>.<div class="im"><br>
Visit this group at <a href="http://groups.google.com/group/cp2k" target="_blank">http://groups.google.com/<u></u>group/cp2k</a>.<br>
For more options, visit <a href="https://groups.google.com/groups/opt_out" target="_blank">https://groups.google.com/<u></u>groups/opt_out</a>.<br>
<br>
<br>
</div></div></div></blockquote></div><br></div>
</blockquote></div><div class=""><div class="h5">
<p></p>
-- <br>
You received this message because you are subscribed to the Google Groups "cp2k" group.<br>
To unsubscribe from this group and stop receiving emails from it, send an email to <a href="mailto:cp2k%2Bun...@googlegroups.com" target="_blank">cp2k+uns...@googlegroups.com</a>.<br>
To post to this group, send email to <a href="mailto:cp...@googlegroups.com" target="_blank">cp...@googlegroups.com</a>.<br>
Visit this group at <a href="http://groups.google.com/group/cp2k" target="_blank">http://groups.google.com/group/cp2k</a>.<br>
For more options, visit <a href="https://groups.google.com/groups/opt_out" target="_blank">https://groups.google.com/groups/opt_out</a>.<br>
<br>
<br>
</div></div></blockquote></div><br></div></div></div>