Re: question about setting up FEP calculations

From: Shyno Mathew (sm3334_at_columbia.edu)
Date: Thu Nov 14 2013 - 15:25:45 CST

Hey Jerome,
the problem was with waters, I wasn't using the equilibrated waters.
I am able to do the FEP calculations.
I have 2 other questions:
1. For doing the backward FEP calculations ( i.e lambda 1 to 0), I can use
the same starting configurations as used for forward FEP, correct?
Let's say I am mutating A --> B
Ideally the backward calculations should give me the same results as the
forward FEP (ie. lambda 0 to 1) but with different starting configuration?

2. Is there a rule of thumb for choosing delta lambda and length of FEP
calculations?
I am using dlamba=0.05
For each lambda window, I used 5000 steps minimization followed by 10,000
step run.

thanks,
Shyno

On Sat, Nov 2, 2013 at 5:24 AM, Jérôme Hénin <jerome.henin_at_ibpc.fr> wrote:

> Well, something is off with either your starting coordinates or the
> topology. It's hard to say anything more precise based on what you said.
>
> What I'd do now is run a long minimization (say at least 10K steps) and
> inspect the minimization trajectory carefully in VMD.
>
> Cheers,
> Jerome
>
>
> ----- Original Message -----
> > Hey Jerome,
> > I realized the co-ordinates of the appearing atoms I was using (for the
> > previous discussions) were very far off. So I wrote a new script and made
> > the initial files which has all atoms (appearing, disappearing and the
> > unchanging ones).
> > Again, I am getting similar errors as before. Also, I did use alchemify
> on
> > the initial files.
> > I changed the time step from 1 fs to 0.003 fs. Also tried playing with
> the
> > cutoff, pairlist dist.
> > Finally time step =0.003 fs seems to work (before it wasn't passing the
> > minimization stage)
> > After 3100 steps, I am getting the error:
> > FATAL ERROR: Low global exclusion count! System unstable or pairlistdist
> > or cutoff too small.
> >
> > I am attaching the output file, if that helps.
> >
> > Following are the info in my .conf file:
> > readexclusions no
> > exclude scaled1-4
> > 1-4scaling 0.833333333
> > scnb 2.0
> > temperature $temperature
> >
> > # OUTPUT FREQUENCIES
> >
> > outputenergies 100
> > outputtiming 100
> > outputpressure 100
> > restartfreq 100
> > XSTFreq 100
> >
> >
> > # OUTPUT AND RESTART
> >
> > outputname 06collab
> > restartname 06collab
> >
> > binaryoutput yes
> > binaryrestart yes
> >
> > # CONSTANT-T
> > langevin on
> > langevinTemp $temperature
> > langevinDamping 5.0
> >
> > langevinHydrogen off
> >
> >
> >
> > # Periodic Boundary Conditions
> > cellBasisVector1 124.16300010681152 0. 0.
> > cellBasisVector2 0. 82.71599864959717 0.
> > cellBasisVector3 0. 0. 128.3490011692047
> > cellOrigin -20.06950092315674 -22.470999240875244
> 46.29550063610077
> >
> >
> > # PME
> > PME yes
> > PMEGridSizeX 128
> > PMEGridSizeY 96
> > PMEGridSizeZ 160
> >
> > # WRAP WATER FOR OUTPUT
> >
> > wrapAll on
> >
> > # CONSTANT-P
> >
> > langevinPiston on
> > langevinPistonTarget 1.01325 ;# in bar -> 1 atm
> > langevinPistonPeriod 100.
> > langevinPistonDecay 50.
> > langevinPistonTemp $temperature
> >
> >
> > #StrainRate 0.0 0.0 0.0
> > useGroupPressure yes
> >
> > useFlexibleCell no
> > useConstantArea no
> >
> > # SPACE PARTITIONING
> > stepspercycle 10
> > margin 3.0
> >
> > # CUT-OFFS
> >
> > switching on
> > switchdist 10.0
> > cutoff 12.0
> > pairlistdist 13.5
> >
> >
> > # RESPA PROPAGATOR
> >
> > timestep 0.003
> > fullElectFrequency 2
> > nonbondedFreq 1
> >
> >
> > # SHAKE
> >
> > rigidbonds all
> >
> >
> >
> > # FEP PARAMETERS
> >
> >
> > alch on
> > alchType FEP
> > alchFile 05collab_wet.fep
> > alchCol B
> > alchOutFile 06collab.fepout
> > alchOutFreq 100
> > alchVdwLambdaEnd 1.0
> > alchElecLambdaStart 0.5
> > alchVdWShiftCoeff 5.0
> > alchDecouple off
> > alchEquilSteps 2000
> >
> > set Lambda0 0.0
> > set dLambda 0.05
> > while {$Lambda0 < 1.0} {
> > alchLambda $Lambda0
> > set Lambda0 [expr \$Lambda0 + \$dLambda]
> > alchLambda2 $Lambda0
> > run 2500
> > }
> >
> >
> >
> > thanks again,
> > Shyno
> >
> >
> >
> > On Fri, Oct 18, 2013 at 12:44 PM, Shyno Mathew <sm3334_at_columbia.edu>
> wrote:
> >
> > > Thanks again Jerome. I did try time steps 0.5, 0.25 still the same
> error
> > > persists.
> > > Please bear with me as I explain what I did:
> > >
> > > 1. Changed time step to 0.5
> > > Reason: FATAL ERROR: Low global exclusion count! System unstable or
> > > pairlistdist or cutoff too small.
> > > FATAL ERROR: See http://www.ks.uiuc.edu/Research/namd/bugreport.html
> > > Charm++ fatal error:
> > > FATAL ERROR: Low global exclusion count! System unstable or
> pairlistdist
> > > or cutoff too small.
> > > FATAL ERROR: See http://www.ks.uiuc.edu/Research/namd/bugreport.html
> > >
> > > 2. So I changed the pairlistdist, cutoff as follows:
> > > time step =0.5, margin =3
> > > switchdist 15.0
> > > cutoff 16.0
> > > pairlistdist 17.5
> > > Reason: FATAL ERROR: Periodic cell has become too small for original
> patch
> > > grid!
> > > Possible solutions are to restart from a recent checkpoint,
> > > increase margin, or disable useFlexibleCell for liquid simulation.
> > >
> > > 3. increased margin upto 6, not sure if this makes sense
> > > Then I increased margin, but still there is an error
> > > Charm++ fatal error:
> > > FATAL ERROR: Periodic cell has become too small for original patch
> grid!
> > > Possible solutions are to restart from a recent checkpoint,
> > > increase margin, or disable useFlexibleCell for liquid simulation.
> > >
> > > 4. Description of the system: Disappearing atom is Hydrogen attached
> to a
> > > carbon atom. Then a methyl group (connected to the previous carbon
> atom )
> > > appears. In the initial file where I have all the appearing and
> > > disappearing atoms, the appearing atoms (methyl group) apppears to be
> very
> > > far from the carbon atom it is connected. I am thinking of minimizing
> this
> > > initial set up first and then submitting the FEP calculation?
> > >
> > > thanks,
> > > Shyno
> > >
> > >
> > >
> > > On Fri, Oct 18, 2013 at 5:33 AM, Jérôme Hénin <jerome.henin_at_ibpc.fr
> >wrote:
> > >
> > >> What happens if you start the MD with a 0.5 fs time step?
> > >>
> > >> Jerome
> > >>
> > >> ----- Original Message -----
> > >> > Thanks again Jerome.
> > >> > As you said I ran alchemify first on the psf file. Here is the
> output
> > >> info
> > >> > I got:
> > >> >
> > >> > ~/Alchemify/alchemify_LINUX 05collab_wet.psf out.psf
> 05collab_wet.fep
> > >> >
> > >> > FEPfile : 112454 atoms found, 10 initial, 40 final
> > >> > 69937 angles : removing 0 angles coupling initial and final groups
> > >> > 57708 dihedrals : removing 0 dihedrals coupling initial and final
> groups
> > >> > 3932 impropers : removing 0 impropers coupling initial and final
> groups
> > >> > Writing 400 exclusion pairs
> > >> >
> > >> > Then I tried running FEP calculations. But I am still getting the
> > >> following
> > >> > error:
> > >> >
> > >> > Reason: FATAL ERROR: Low global exclusion count! System unstable or
> > >> > pairlistdist or cutoff too small.
> > >> >
> > >> > FATAL ERROR: See
> http://www.ks.uiuc.edu/Research/namd/bugreport.html
> > >> >
> > >> > Charm++ fatal error:
> > >> > FATAL ERROR: Low global exclusion count! System unstable or
> > >> pairlistdist
> > >> > or cutoff too small.
> > >> >
> > >> >
> > >> > I tried playing with cutoff distance!
> > >> >
> > >> > thanks again for your help,
> > >> > Shyno
> > >> >
> > >> >
> > >> >
> > >> >
> > >> > On Tue, Oct 15, 2013 at 1:00 PM, Jérôme Hénin <jerome.henin_at_ibpc.fr
> >
> > >> wrote:
> > >> >
> > >> > > No, that short run is not correct, because it does not describe
> any
> > >> > > physically meaningful state.
> > >> > >
> > >> > > Now I suspect one thing: that the two end-point groups are
> clashing,
> > >> and
> > >> > > that causes instability. The only reason why that could happen is
> if
> > >> you
> > >> > > don't have the proper nonbonded exclusions, which could happen if
> you
> > >> are
> > >> > > running the memory-optimized version of NAMD.
> > >> > >
> > >> > > To make sure, you could try running Alchemify on your PSF before
> > >> running:
> > >> > > http://www.edam.uhp-nancy.fr/Alchemify/
> > >> > >
> > >> > > Best,
> > >> > > Jerome
> > >> > >
> > >> > >
> > >> > > ----- Original Message -----
> > >> > > > thanks again Jerome. Sorry my question was kind of hidden .
> > >> > > > Just to be clear,
> > >> > > > If I submit FEP calculations straight from the dual topology
> file,
> > >> it
> > >> > > > wasn't running. So when I submitted the short run (~ 3ns) and
> then
> > >> did
> > >> > > FEP
> > >> > > > calculations it worked.
> > >> > > > So, as mentioned before, for this short run, the system has both
> > >> hydrogen
> > >> > > > (disappearing atom) and methyl group (appearing) bonded to the
> > >> carbon
> > >> > > atom
> > >> > > > at the same time. I wasn't sure if this part (short run with all
> > >> > > > atoms(appearing & disappearing ) before FEP ) is correct or not?
> > >> > > >
> > >> > > > thanks,
> > >> > > > Shyno
> > >> > > >
> > >> > > >
> > >> > > > On Tue, Oct 15, 2013 at 12:29 PM, Jérôme Hénin <
> > >> jerome.henin_at_ibpc.fr>
> > >> > > wrote:
> > >> > > >
> > >> > > > > Sorry, there was a misunderstanding. Then what you were doing
> > >> seems
> > >> > > > > correct, and I don't know why it didn't run properly.
> > >> > > > >
> > >> > > > > Jerome
> > >> > > > >
> > >> > > > >
> > >> > > > > ----- Original Message -----
> > >> > > > > > thanks again for your reply. Sorry for asking more
> questions.
> > >> Please
> > >> > > bear
> > >> > > > > > with me as I explain few details:
> > >> > > > > > This is exactly what I did. The dual topology was made in
> the
> > >> > > beginning
> > >> > > > > (i
> > >> > > > > > used non equilibrated systems to make it). Then added
> > >> counterions and
> > >> > > > > > water. After this the dual topology system was minimized and
> > >> > > submitted
> > >> > > > > for
> > >> > > > > > a short run before FEP calculations. As mentioned before,
> for
> > >> this
> > >> > > short
> > >> > > > > > run, the system has both hydrogen (disappearing atom) and
> methyl
> > >> > > group
> > >> > > > > > (appearing) bonded to the carbon atom at the same time
> during
> > >> the 3
> > >> > > ns
> > >> > > > > run.
> > >> > > > > >
> > >> > > > > > Or you meant to say the dual topology has to be made using
> the
> > >> > > > > equilibrated
> > >> > > > > > systems?
> > >> > > > > >
> > >> > > > > > thanks again,
> > >> > > > > > Shyno
> > >> > > > > >
> > >> > > > > >
> > >> > > > > >
> > >> > > > > > On Tue, Oct 15, 2013 at 11:44 AM, Jérôme Hénin <
> > >> jerome.henin_at_ibpc.fr
> > >> > > >
> > >> > > > > wrote:
> > >> > > > > >
> > >> > > > > > > In general, it is simpler to setup the complete "dual
> > >> topology"
> > >> > > system
> > >> > > > > > > from the beginning, then add solvent and other components.
> > >> > > > > > >
> > >> > > > > > > Jerome
> > >> > > > > > >
> > >> > > > > > >
> > >> > > > > > > ----- Original Message -----
> > >> > > > > > > > Hello Jerome,
> > >> > > > > > > > Thank you very much for your reply.
> > >> > > > > > > > Ok, so I should use the equilibrated system for lambda=0
> > >> (this is
> > >> > > > > where I
> > >> > > > > > > > want to start). To this I can add the appearing atoms
> > >> (lambda=1)
> > >> > > also
> > >> > > > > > > after
> > >> > > > > > > > the equilibration run and then create the .fep file and
> so
> > >> on.
> > >> > > > > > > > Initially, I was trying this approach, one difficulty I
> > >> faced was
> > >> > > > > how to
> > >> > > > > > > > number the atoms after 99999.
> > >> > > > > > > > For example, when I add the appearing atoms to the pdb
> file
> > >> > > > > (lambda=0 )
> > >> > > > > > > > after equilibration, the seriel numbers change. This
> problem
> > >> > > arises
> > >> > > > > from
> > >> > > > > > > > the water molecules. Maybe another option would be just
> to
> > >> select
> > >> > > > > > > > everything except water (for lambda=0 after
> equilibration)
> > >> and
> > >> > > add
> > >> > > > > > > > appearing atoms and then solvate the system? Is this
> > >> approach
> > >> > > > > correct?
> > >> > > > > > > >
> > >> > > > > > > > thanks again for your help,
> > >> > > > > > > > Shyno
> > >> > > > > > > >
> > >> > > > > > > >
> > >> > > > > > > >
> > >> > > > > > > > On Mon, Oct 14, 2013 at 3:32 PM, Jérôme Hénin <
> > >> > > jerome.henin_at_ibpc.fr>
> > >> > > > > > > wrote:
> > >> > > > > > > >
> > >> > > > > > > > > Hi Shyno,
> > >> > > > > > > > >
> > >> > > > > > > > > There seems to be something wrong with your approach--001a11c345f687328f04eb29b912--

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2013 - 23:23:59 CST