I mean the q_tot^2 term. This attractative term goes with the square of the total charge, so has larger magnitude for AB than A+B, so the result in negative.
Berk ________________________________ From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se <gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Sergio Perez <sperezcon...@gmail.com> Sent: Friday, March 29, 2019 3:15 PM To: gmx-us...@gromacs.org Subject: Re: [gmx-users] PME constant shifts. Hello, Thanks for the quick reply, I am not sure if you are referring to this term of the Ewald Sum: V_0 = ( f*beta ) * SUM(q_i^2) / sqrt(pi) Since the sum of the square of the charges of the system (AB) is the same as the sum of the square of the charges of the system (A) plus the same from system (B), this term should vanish. There is this an additional term for the Ewald Sum when the system is charged to compensate for the uniform charge as described by K. Smith in Elements of Molecular Dynamics section 6.4.6 (free book online): V= - (q_tot^2)/(8*eps_0*V*Beta) If this was not included by gromacs then this would explain why I get negative interaction energies. I know charged systems should be avoided (for example if the system is inhomogeneous, as grompp points out) but my purposes are FF development and finding bugs in my FF. Although the absolute energies are meaningless it would be nice if they could be intuitively correct. Best regards and thanks again! Sergio On Fri, Mar 29, 2019 at 1:45 PM Berk Hess <g...@hotmail.com> wrote: > Hi, > > The negative energy is due to the uniform background charge that appears > when your system has a net charge, not due to shifting of pair > interactions. Computing the absolute energy of a periodic system with net > charge is meaningless. > > Cheers, > > Berk > > ________________________________ > From: gromacs.org_gmx-users-boun...@maillist.sys.kth.se < > gromacs.org_gmx-users-boun...@maillist.sys.kth.se> on behalf of Sergio > Perez <sperezcon...@gmail.com> > Sent: Friday, March 29, 2019 12:33 PM > To: gromacs.org_gmx-users@maillist.sys.kth.se > Subject: [gmx-users] PME constant shifts. > > Dear gmx comunity, > > I am trying to calculate the electrostatic interaction of my system for > force field development. My system consists of 7 positively charged > particles interacting (system A) with another positively charged particle > (system B). I calculate the electrostatic interaction by doing reruns of > the trajectory in which some of them have been removed. > > E_int(electrostatic) = > > E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR) > > And to my surprise I get negative energies. I know the PME method has its > energy shifted: > > "With the Verlet cut-off scheme, the PME direct space potential is shifted > by a constant such that > the potential is zero at the cut-off. This shift is small and since the net > system charge is close to > zero, the total shift is very small, unlike in the case of the > Lennard-Jones potential where all shifts > add up. We apply the shift anyhow, such that the potential is the exact > integral of the force." > > Since for the system with just one charge I get a possitive E_1q(coul-LR), > I imagine the LR-coul term is shifted (eventhough I have to use the group > cut-off scheme and not the verlet). The problem is why don't all these > shifts cancel? Is there a way to not add the shifts? > > I leave bellow the .mdp parameters. Note that for the VDW part I need to > use user-tables and therefore I need to use the group scheme. > > Thank you very much! > Sergio PĂ©rez-Conesa > > > ; Integrator stuff > integrator = md-vv ; steep or md > dt = 0.001 > nsteps = 300000000 ; -1 no maximum > init-step = 30305000 ; > ld-seed = -2018962629 > continuation = yes > emtol = 10.0 > emstep = 0.01 > > ; neighbour search > cutoff-scheme = group > nst-list = 10 > verlet-buffer-tolerance = 0.005 > ns-type = grid > ; Pbc > pbc = xyz ; xyz or xy > > ; Vdw > rvdw = 1.4 > rlist = 1.4 > vdwtype = user ; PME or User to look for table /user > DispCorr = No ; corrections to energy and pressure or No > > ;Electrostatic > coulombtype = PME ; User if look for table > rcoulomb = 1.4 > pme-order = 4 > fourierspacing = 0.2 > ewald-geometry = 3d ; 3d or 3dc > > ;IMD-group = > > ; T and P > tcoupl =no ; none, nose-hoover, v-rescale > nh-chain-length = 10 > nsttcouple = 10 ; 10 if not used > pcoupl = no ; > pcoupltype = semiisotropic ; > compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility of > water, bar^-1 > refcoord-scaling = com > tau-p = 5.0 ; ps > ref-p = 1.0 1.0 ; bar > > ; Constraints > constraints = none ; constraints distances > constraint-algorithm = LINCS ; or SHAKE > lincs_iter = 1 > lincs_order = 4 > > ;Generate velocities > gen_vel = no ; assign velocities from Maxwell > distribution > gen_temp = 300 ; temperature for Maxwell distribution > gen_seed = 180611 ; generate a random seed > > ; OUTPUT > ; Output control > nstenergy = 1 ; save energies > nstlog = 10 ; update log file every > ;nstxout-compressed = 500 ; save compressed coordinates every 10.0 ps > ;nstxout = 20 > ;nstvout = 20 > ;nstfout = 20 > > ; COM removal . There are more options like removing from groups or the > frequency > ;comm-grps = System > comm-mode = none ; linear, angular (removes linear and angular motion) > or none > nstcomm = 0 ; frequency of removal of com motion > freezegrps = > energygrp-excl = > ;freezedim = Y Y Y > nstcalcenergy = > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.