Re: question about setting up FEP calculations

From: Shyno Mathew (sm3334_at_columbia.edu)
Date: Fri Nov 01 2013 - 10:14:36 CDT

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. If
>> you
>> > > do an
>> > > > > > > > > equilibration run prior to the FEP, you should do that at
>> a
>> > > > > > > well-defined
>> > > > > > > > > lambda state, by enabling the alchemical options and
>> setting
>> > > > > lambada to
>> > > > > > > > > either 0 or 1, depending on where you want to start from--001a11c3145c40c06a04ea1f067f--
--001a11c3145c40c06d04ea1f0681--

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:51 CST