[gmx-users] PME constant shifts.

Sergio Perez sperezconesa at gmail.com
Fri Mar 29 16:25:05 CET 2019 

Great! So I just have to discount this term to obtain the positive energy
expected.
Thank you very much!
Sergio

On Fri, Mar 29, 2019 at 4:11 PM Berk Hess <gmx3 at hotmail.com> wrote:

> 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-bounces at maillist.sys.kth.se <
> gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Sergio
> Perez <sperezconesa at gmail.com>
> Sent: Friday, March 29, 2019 3:15 PM
> To: gmx-users at 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 <gmx3 at 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-bounces at maillist.sys.kth.se <
> > gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Sergio
> > Perez <sperezconesa at gmail.com>
> > Sent: Friday, March 29, 2019 12:33 PM
> > To: gromacs.org_gmx-users at 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-request at 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-request at 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-request at 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-request at gromacs.org.
>

More information about the gromacs.org_gmx-users mailing list