Re: [AMBER] Problems with ti using DFTB for zinc bound hydroxyl

From: Brian Radak <radak004.umn.edu>
Date: Mon, 4 Feb 2013 14:51:45 -0500

What do your input files and group file look like? Also, if you are using
the normal TI code you must provide prmtop files with identical numbers of
atoms.

On Mon, Feb 4, 2013 at 2:21 PM, Fabrício Bracht <bracht.iq.ufrj.br> wrote:

> Ok. I removed the softcore potential and tried to run the simulation for
> lambda 1.0 but the simulation crashes before it even starts. The error
>
> Running multisander version of sander Amber12
> Total processors = 2
> Number of groups = 2
>
> [woodstock:18042] *** An error occurred in MPI_Sendrecv
> [woodstock:18042] *** on communicator MPI COMMUNICATOR 4 SPLIT FROM 0
> [woodstock:18042] *** MPI_ERR_TRUNCATE: message truncated
> [woodstock:18042] *** MPI_ERRORS_ARE_FATAL (your MPI job will now abort)
> --------------------------------------------------------------------------
> mpirun has exited due to process rank 1 with PID 18042 on
> node woodstock exiting without calling "finalize". This may
> have caused other processes in the application to be
> terminated by signals sent by mpirun (as reported here).
>
> Comes out on screen. The command given was
> mpirun -np 2 sander.MPI -ng 2 -groupfile group_heat1
>
> Am I doing something wrong here?
>
> 2013/2/4 Brian Radak <radak004.umn.edu>
>
> > Fabricio,
> >
> > Unfortunately, I have no idea if the softcore potential code is
> compatible
> > with the QM/MM code. I believe the topology information is handled in a
> > slightly unusual way, which may cause trouble. Perhaps Thomas (or Ross?)
> > will see this and be able to comment. However I can't imagine it is
> > necessary to use the softcore code since the proton likely has no
> > Lennard-Jones term anyway. Even if that is not the case, I would guess
> that
> > the epsilon parameters are probably small enough that they can be scaled
> > away linearly, although this is pure speculation.
> >
> > I would try using the normal TI code and setting lambda = 1 or 0 and
> > verifying that you obtain the same dynamics you did with the QM/MM code
> on
> > its own.
> >
> > Regards,
> > Brian
> >
> > On Mon, Feb 4, 2013 at 1:31 PM, Fabrício Bracht <bracht.iq.ufrj.br>
> wrote:
> >
> > > What I meant is that I have already performed a MM ti simulation, but
> > since
> > > I am using soft core potentials, I cannot set Lambda = 1 or lambda = 0.
> > > This means that I have not done such simulations (lambda = 0 or 1) with
> > > qmmm either. Setting restraints to the bonds between the zinc atom and
> > the
> > > hydroxyl group has not worked either. The addition of these restraints
> > > result in very abnormal behaviour, such as very distorted geometry of
> the
> > > entire qm region, as well as SCF failing to converge. What I still
> don't
> > > understand is why such abnormal behaviour happens only with TI
> > simulations.
> > > As I said before, if I simulate only the hydroxyl system (no
> multisander
> > > simulation with the HOH system) with qmmm, the simulation runs
> smoothly.
> > Do
> > > you imagine I will have more success doing Ti simulation without
> softcore
> > > potentials?
> > > Fabrício Bracht
> > >
> > > 2013/2/3 Brian Radak <radak004.umn.edu>
> > >
> > > > Fabricio,
> > > >
> > > > So do you mean that setting lambda = 1.0 or 0.0 with the TI code
> gives
> > > > different results than the QM/MM code on its own? That would be
> > > concerning
> > > > indeed. I'm pretty sure that, all else equal, running an NVE
> simulation
> > > > (ntt = 0) or setting ig to an explicit value should result in
> identical
> > > > trajectories for both of those circumstances, even with multiple
> > > > processors. That is the most simplest and most obvious test that I
> can
> > > > think of (nstlim = 10, ntpr = 1, etc. should be adequate).
> > > >
> > > > As for the alternate transformation method, I am not sure how much
> this
> > > is
> > > > done in the literature, but the idea would be to have the proton
> > > > "disappear" via an MM transformation with an appropriate force field
> > (in
> > > > principle this can be arbitrary). Then, one would perform two
> > additional
> > > > transformations from the MM endpoints where one endpoint is the MM
> > > > description and the other is the analogous QM/MM description. This
> > should
> > > > in principle amount to a smaller perturbation of the bonded and
> > > > electrostatic terms (no guarantees!) as well as avoid the "problem"
> of
> > > > having weird bonded terms on the QM atoms.
> > > >
> > > > I have no idea if the latter suggestion would solve your instability
> > > > issues, but if no other solutions present themselves it's probably
> > worth
> > > a
> > > > try.
> > > >
> > > > Regards,
> > > > Brian
> > > >
> > > > On Sat, Feb 2, 2013 at 8:33 PM, Fabrício Bracht <bracht.iq.ufrj.br>
> > > wrote:
> > > >
> > > > > Hi Brian. Thank you for the advice. Indeed I had already seen the
> > > > article,
> > > > > but I think it is time to read it again with more attention. The
> > lambda
> > > > > values I am using are 0.01592; 0.08198; 0.19331; 0.33787; 0.5;
> > 0.66213;
> > > > > 0.80669; 0.91802; 0.98408. These are the values listed in amber12
> > > manual
> > > > > page #117.
> > > > > The odd thing though is the difference I get in behaviour when I
> run
> > > the
> > > > > simulation still with qmmm and using the same basic input
> parameters;
> > > > i.e.,
> > > > > using the same input file as for the TI calculation, but without
> the
> > TI
> > > > > part. All atoms stay bound, no matter what happens. At least for
> the
> > > > > hydroxyl ligand. I was, indeed, expecting some problems to arise
> with
> > > the
> > > > > water molecule, since an oxygen atom bound to two hydrogens and a
> > zinc
> > > > atom
> > > > > is unlikely to stay that way for long.
> > > > > You mentioned doing a MM simulation. Well, I already did that,
> but I
> > > am
> > > > > not sure I understood what you meant by "perform MM to QM/MM
> > > > > transformations on the endpoints as a kind of correction."
> > > > > Thank you again
> > > > > Fabrício Bracht
> > > > >
> > > > >
> > > > > 2013/2/2 Brian Radak <radak004.umn.edu>
> > > > >
> > > > > > Hi Fabricio,
> > > > > >
> > > > > > Can you clarify what exactly you mean by "start with the TI
> > > > simulation"?
> > > > > > What values of lambda? Does this behavior occur for all of them?
> > > > > >
> > > > > > Some other things that might be informative: What topolog(y/ies)
> > are
> > > > you
> > > > > > using? What bonded terms exist on the protons?
> > > > > >
> > > > > > I believe the best established QM/MM protocols for pKa type
> > > > calculations
> > > > > > involve some kind of restraints, but how you choose to apply them
> > is
> > > > > still
> > > > > > something of an art. The following paper describes a protocol in
> > > > CHARMM,
> > > > > > but I'm pretty sure the ideas transfer well to AMBER:
> > > > > >
> > > > > > Li and Qiang, J Phys Chem B 2003, 107, 14521.
> > > > > >
> > > > > > I attempted similar calculations a while back but had a very
> > > difficult
> > > > > time
> > > > > > getting reproducible results. I suspect that there are some zero
> of
> > > > > energy
> > > > > > issues that need to be handled when you have the same number of
> > > > electrons
> > > > > > but different orbitals. I am not sure what the solution is there.
> > An
> > > > > > alternative path is to perform the proton disappearing as an MM
> to
> > MM
> > > > > > mutation and then perform MM to QM/MM transformations on the
> > > endpoints
> > > > > as a
> > > > > > kind of correction. This is more work but I suspect will converge
> > > > faster
> > > > > in
> > > > > > most situations.
> > > > > >
> > > > > > Regards,
> > > > > > Brian
> > > > > >
> > > > > > On Sat, Feb 2, 2013 at 2:35 PM, Fabrício Bracht <
> bracht.iq.ufrj.br
> > >
> > > > > wrote:
> > > > > >
> > > > > > > Hi all. I am doing a TI caltulation in order to obtain the free
> > > > energy
> > > > > > > values associated with the loss of a proton from a hydroxyl
> group
> > > > bound
> > > > > > to
> > > > > > > a zinc atom coordinated with 3 other aminoacids. System 1
> > comprises
> > > > of
> > > > > a
> > > > > > > water molecule bound to the zinc atom, and system 2 is the
> > hydroxyl
> > > > > > group.
> > > > > > > When I run qmmm simulations, using DFTB, of both systems
> > > separately,
> > > > > the
> > > > > > > system stays intact, i.e. no bonds are broken and both the
> water
> > > and
> > > > > the
> > > > > > > hydroxyl stay bound to the zinc atom ( as was expected). But as
> > > soon
> > > > > as I
> > > > > > > star with the TI simulation, the hydroxyl group always drifts
> > away.
> > > > > > > Sometimes the oxygen atom stays and the hydrogen goes for a
> walk.
> > > > This
> > > > > > > happens only on TI simulations and the water system stays
> > intact. I
> > > > > have
> > > > > > > tried different sizes of qm region, different cutoffs, tried
> > > several
> > > > > > > different starting structures ( minimized, equilibrated,
> > > thermalized
> > > > > etc)
> > > > > > > and the result is almost always the same. Below is the input
> file
> > > for
> > > > > the
> > > > > > > hydroxyl group system. I thought about restraining the hydroxyl
> > > group
> > > > > > using
> > > > > > > the nmropt module. Can I use the dvdl_norest flag to ignore
> those
> > > > > > > restraints in the energy calculation? Please let me know if you
> > > need
> > > > > any
> > > > > > > other detail from the simulation and/or if you need/want any
> file
> > > > from
> > > > > > the
> > > > > > > inputs.
> > > > > > > Thank you
> > > > > > > Fabrício Bracht
> > > > > > >
> > > > > > > TI of hydroxyl group step 1
> > > > > > > &cntrl
> > > > > > > imin = 0,
> > > > > > > irest = 1,
> > > > > > > ntx = 7,
> > > > > > > ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0,
> > > > > > > cut = 8.0,
> > > > > > > ntr = 0,
> > > > > > > ntc = 2,
> > > > > > > ntf = 1,
> > > > > > > temp0 = 298.0,
> > > > > > > ntt = 3,
> > > > > > > gamma_ln = 1.0,
> > > > > > > nstlim = 1000000, dt = 0.0005, ntave = 100,
> > > > > > > ntpr = 100, ntwx = 100, ntwr = 100,
> > > > > > > ig = 10703, ioutfm = 1, iwrap = 1,
> > > > > > > icfe = 1, ifsc = 1, clambda = 0.01592,
> > > > > > > ifqnt = 1, scmask = ':342', idecomp = 0, nmropt = 0,
> > > > > > > /
> > > > > > > &qmmm
> > > > > > > qmmask=
> > > > > > >
> > > ':200,255,258&!.CA,C,HA,O,N,HN,H|:341,342|(:202,204&!.CA,C,HA,O,N,H)'
> > > > > > > dftb_3rd_order = 'PA'
> > > > > > > qmcharge=-1,
> > > > > > > qm_theory='DFTB',
> > > > > > > qmshake=0,
> > > > > > > qm_ewald=1, qm_pme=1
> > > > > > > qmcut=9.0
> > > > > > > writepdb=1
> > > > > > > /
> > > > > > > Receptor residues
> > > > > > > RRES 1 7853
> > > > > > > END
> > > > > > > Printing
> > > > > > > RES 1 342
> > > > > > > END
> > > > > > > END
> > > > > > > _______________________________________________
> > > > > > > AMBER mailing list
> > > > > > > AMBER.ambermd.org
> > > > > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > --
> > > > > > ================================ Current Address
> > > > =======================
> > > > > > Brian Radak :
> > > BioMaPS
> > > > > > Institute for Quantitative Biology
> > > > > > PhD candidate - York Research Group : Rutgers, The
> State
> > > > > > University of New Jersey
> > > > > > University of Minnesota - Twin Cities : Center for
> > > > > Integrative
> > > > > > Proteomics Room 308
> > > > > > Graduate Program in Chemical Physics : 174 Frelinghuysen
> > > Road,
> > > > > > Department of Chemistry :
> Piscataway,
> > > NJ
> > > > > > 08854-8066
> > > > > > radak004.umn.edu :
> > > > > > radakb.biomaps.rutgers.edu
> > > > > >
> > ====================================================================
> > > > > > Sorry for the multiple e-mail addresses, just use the institute
> > > > > appropriate
> > > > > > address.
> > > > > > _______________________________________________
> > > > > > AMBER mailing list
> > > > > > AMBER.ambermd.org
> > > > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > > > >
> > > > > _______________________________________________
> > > > > AMBER mailing list
> > > > > AMBER.ambermd.org
> > > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > > >
> > > >
> > > >
> > > >
> > > > --
> > > > ================================ Current Address
> > =======================
> > > > Brian Radak :
> BioMaPS
> > > > Institute for Quantitative Biology
> > > > PhD candidate - York Research Group : Rutgers, The State
> > > > University of New Jersey
> > > > University of Minnesota - Twin Cities : Center for
> > > Integrative
> > > > Proteomics Room 308
> > > > Graduate Program in Chemical Physics : 174 Frelinghuysen
> Road,
> > > > Department of Chemistry : Piscataway,
> NJ
> > > > 08854-8066
> > > > radak004.umn.edu :
> > > > radakb.biomaps.rutgers.edu
> > > > ====================================================================
> > > > Sorry for the multiple e-mail addresses, just use the institute
> > > appropriate
> > > > address.
> > > > _______________________________________________
> > > > AMBER mailing list
> > > > AMBER.ambermd.org
> > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > >
> > > _______________________________________________
> > > AMBER mailing list
> > > AMBER.ambermd.org
> > > http://lists.ambermd.org/mailman/listinfo/amber
> > >
> >
> >
> >
> > --
> > ================================ Current Address =======================
> > Brian Radak : BioMaPS
> > Institute for Quantitative Biology
> > PhD candidate - York Research Group : Rutgers, The State
> > University of New Jersey
> > University of Minnesota - Twin Cities : Center for
> Integrative
> > Proteomics Room 308
> > Graduate Program in Chemical Physics : 174 Frelinghuysen Road,
> > Department of Chemistry : Piscataway, NJ
> > 08854-8066
> > radak004.umn.edu :
> > radakb.biomaps.rutgers.edu
> > ====================================================================
> > Sorry for the multiple e-mail addresses, just use the institute
> appropriate
> > address.
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>



-- 
================================ Current Address =======================
 Brian Radak                                             :     BioMaPS
Institute for Quantitative Biology
 PhD candidate - York Research Group       :     Rutgers, The State
University of New Jersey
 University of Minnesota - Twin Cities         :     Center for Integrative
Proteomics Room 308
 Graduate Program in Chemical Physics     :     174 Frelinghuysen Road,
 Department of Chemistry                          :     Piscataway, NJ
08854-8066
 radak004.umn.edu                                 :
radakb.biomaps.rutgers.edu
====================================================================
Sorry for the multiple e-mail addresses, just use the institute appropriate
address.
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Feb 04 2013 - 12:00:03 PST
Custom Search