[gmx-users] Free Energy solvation generate multiple dhdl.xvg output files

Javier Luque Di Salvo jluquedisalvo at gmail.com
Tue Mar 20 13:18:19 CET 2018


Dear Wes,

Thanks for clarifying, I did as you told and it worked. The resulting
dG_hyd of NH4+ is -62.5 kJ/mol, while the experimental value is -285 kJ/mol
(10.1039/FT9918702995). I am using OPLS and SPC/E water. Energy values of
the last lambdas showed higher error bars, so now I am refining the last
part by increasing the number of lambda points:
fep-lambdas              = 0.0 0.2 0.4 0.6 0.8 0.85 0.9 0.92 0.94 0.96 0.98
1.0
I do not expect big differences. I could also try (i) other water model,
(ii) bigger box with more NH4+ ions (same concentration) or (ii-b) increase
NH4+ concentration, (iii) adding anions (like Cl-). Regarding these
options, I am not sure yet if they have sense (meaning that there may be
another more suitable option, while some of the possible tests I wrote may
not result in significative improvements), so if there are any comments
regarding these options they are welcome.

Also, I do not understand why the Nitrogen partial charge is negative,
shouldnt be positive? I choose OPLS since it was parameterized for liquid
simulations, but I found in literature that usually there are problems for
correctly simulating dG_hyd of amines.
 atom     charge       mass
 N00    -0.5183    14.0070
 H01     0.3796     1.0080
 H02     0.3796     1.0080
 H03     0.3796     1.0080
 H04     0.3796     1.0080
Maybe the best option is to change FF (does someone know about an
acceptable FF -preferably all_atoms- for amines/ ammonium derivatives R4N+)
?

Best regards
Javier

Date: Mon, 19 Mar 2018 13:50:19 -0400
> From: Wes Barnett <w.barnett at columbia.edu>
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Free Energy solvation generate multiple
>         dhdl.xvg output files
> Message-ID:
>         <CALkN6e3i4ity-kZOf5c7eZO=5jumdrBPz620LQK_1nr2YSExxw at mail.
> gmail.com>
> Content-Type: text/plain; charset="UTF-8"
>
> It sounds like you are only running one free energy simulation. You need to
> run one simulation for each lambda (changing "init-lambda-state" in your
> mdp file) and name the dhdl xvg file a different name for each one.
>
> On Mon, Mar 19, 2018 at 8:30 AM, Javier Luque Di Salvo <
> jluquedisalvo at gmail.com> wrote:
>
> > Dear Gromacs users,
> > I am trying to run Free Energy calculations of a single ion in water (to
> > calculate its solvation energy). Since this is a new topic for me, I am
> > following the tutorial 'Solvation free energy of ethanol' (easily find in
> > google). The three steps work well without erros (1-energy minimization,
> > 2-equilibration, 3-lambda (de)coupling) but I am not obtaining the
> separate
> > dhdl.xvg files after the step 3, only a single output file run.xvg. The
> > reason for this is that in the tutorial, a script named mklambdas.sh is
> > used to generate such separate files in separate folders. Since I would
> > like to keep track of the .mdp strings and understand how they work, I am
> > not using this script, instead I am trying to properly define mdp
> commands:
> >
> > _________________________________________
> > ; Free energy variables
> > free-energy              = yes
> > couple-moltype           = NH4
> > couple-lambda0           = none
> > couple-lambda1           = vdwq
> > couple-intramol          = no
> > init-lambda              = -1
> > init-lambda-state        = 0
> > delta-lambda             = 0
> > nstdhdl                  = 50
> > fep-lambdas              = 0.0 0.2 0.4 0.6 0.8 0.9 1.0
> > mass-lambdas             =
> > coul-lambdas             =
> > vdw-lambdas              =
> > bonded-lambdas           =
> > restraint-lambdas        =
> > temperature-lambdas      =
> > calc-lambda-neighbors    = 1
> > init-lambda-weights      =
> > dhdl-print-energy        = no
> > sc-alpha                 = 1.0
> > sc-power                 = 1
> > sc-r-power               = 6
> > sc-sigma                 = 0.3
> > sc-coul                  = no
> > separate-dhdl-file       = yes
> > dhdl-derivatives         = yes
> > dh_hist_size             = 0
> > dh_hist_spacing          = 0.1
> > _________________________________________
> >
> > When trying to use "gmx bar -f run.xvg -g run.edr -o bar.xvg" I get this
> > error:
> > _______________________________________________________
> > run.xvg: Ignoring set 'pV (kJ/mol)'.
> > run.xvg: 0.0 - 200.0; lambda = 0
> >     dH/dl & foreign lambdas:
> >         dH/dl (fep-lambda) (2001 pts)
> >         delta H to 0 (2001 pts)
> >         delta H to 0.2 (2001 pts)
> >
> > Opened run.edr as single precision energy file
> > Reading energy frame      0 time    0.000
> > -------------------------------------------------------
> > Program:     gmx bar, version 2016.4
> > Source file: src/gromacs/gmxana/gmx_bar.cpp (line 3286)
> > Fatal error:
> > Did not find delta H information in file run.edr
> > ____________________________________________________
> >
> > This is the run.xvg output file:
> >
> > ____________________________________________________
> > # Command line:
> > #   gmx mdrun -ntmpi 2 -ntomp 8 -pin on -v -deffnm run
> > # gmx mdrun is part of G R O M A C S:
> > #
> > # S  C  A  M  O  R  G
> > #
> > @    title "dH/d\xl\f{} and \xD\f{}H"
> > @    xaxis  label "Time (ps)"
> > @    yaxis  label "dH/d\xl\f{} and \xD\f{}H (kJ/mol [\xl\f{}]\S-1\N)"
> > @TYPE xy
> > @ subtitle "T = 300 (K) \xl\f{} state 0: fep-lambda = 0.0000"
> > @ view 0.15, 0.15, 0.75, 0.85
> > @ legend on
> > @ legend box on
> > @ legend loctype view
> > @ legend 0.78, 0.8
> > @ legend length 2
> > @ s0 legend "dH/d\xl\f{} fep-lambda = 0.0000"
> > @ s1 legend "\xD\f{}H \xl\f{} to 0.0000"
> > @ s2 legend "\xD\f{}H \xl\f{} to 0.2000"
> > @ s3 legend "pV (kJ/mol)"
> > 0.0000 -530.14746 0.0000000 -106.84067 9.7838984
> > 0.1000 -177.17302 0.0000000 -33.945797 9.7844772
> > 0.2000 -158.46501 0.0000000 -30.541293 9.7854595
> > 0.3000 -143.72632 0.0000000 -27.669909 9.7857714
> > (...)
> >
> > ____________________________________________________
> >
> >
> > Could you guide me through a proper string definition in the mdp options
> to
> > get separate files, in order to use the BAR algorithm afterwards?
> >
> > Best regards
> >
> > *Javier*
> > --
> > 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.
> >
>
>
>
> --
> James "Wes" Barnett
> Postdoctoral Research Scientist
> Department of Chemical Engineering
> Kumar Research Group <http://www.columbia.edu/cu/kumargroup/>
> Columbia University
> w.barnett at columbia.edu
> http://wbarnett.us
>
>
> ------------------------------
>
> --
> 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.
>
> End of gromacs.org_gmx-users Digest, Vol 167, Issue 93
> ******************************************************
>



-- 

________________________________________

*Javier Luque Di Salvo*

Dipartamento di Ingegneria Chimica

Universitá Degli Studi di Palerm*o*


More information about the gromacs.org_gmx-users mailing list