From: Jackson Cavett (jcavett_at_mail.bradley.edu)
Date: Sat Aug 13 2016 - 11:29:31 CDT

Also, this is the configuration file for the equilibration process. I am
wondering if some of my parameters are wrong within this file. The system
has 170,000 atoms, so it was run for 50 nanoseconds. Please let me know if
there are any glaring issues.

#############################################################
## JOB DESCRIPTION ##
#############################################################

# Equilibration of
# Integrin alphaVbeta3 in a Water Box

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################

structure alphaVbeta3_wb.psf
coordinates alphaVbeta3_wb.pdb

set temperature 310
set outputname PLOS_equil_part_10

# Continuing a job from the restart files
if {1} {
set inputname PLOS_equil_part_9
binCoordinates $inputname.restart.coor
#binVelocities $inputname.restart.vel
extendedSystem $inputname.xsc
}

firsttimestep 190000

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters ./par_all27_prot_lipid.inp
temperature $temperature

# Force-Field Parameters
exclude scaled1-4
1-4scaling 1.0
cutoff 12.0
switching on
switchdist 10.0
pairlistdist 14.0

# Integrator Parameters
timestep 2.0 ;# 2fs/step
rigidBonds all ;# needed for 2fs steps
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 10

# Constant Temperature Control
langevin on ;# do langevin dynamics
langevinDamping 5 ;# damping coefficient (gamma) of 1/ps
langevinTemp $temperature
langevinHydrogen off ;# don't couple langevin bath to hydrogens

# Periodic Boundary Conditions
cellBasisVector1 129.0 0. 0.0
cellBasisVector2 0.0 98.9 0.0
cellBasisVector3 0.0 0 145.3
cellOrigin -16.2 33.9 38.66

wrapAll on

# PME (for full-system periodic electrostatics)
PME yes
PMEGridSpacing 1.0

#manual grid definition
#PMEGridSizeX 45
#PMEGridSizeY 45
#PMEGridSizeZ 48

# Constant Pressure Control (variable volume)
useGroupPressure yes ;# needed for rigidBonds
useFlexibleCell no
useConstantArea no

langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
langevinPistonPeriod 100.0
langevinPistonDecay 50.0
langevinPistonTemp $temperature

# Output
outputName $outputname

restartfreq 500 ;# 500steps = every 1ps
dcdfreq 250
xstFreq 250
outputEnergies 100
outputPressure 100

#############################################################
## EXTRA PARAMETERS ##
#############################################################

#############################################################
## EXECUTION SCRIPT ##
#############################################################

# Minimization
#minimize 100
reinitvels $temperature

run 25000000 ;# 50ns

Thanks once again!

-Jack

On Sat, Aug 13, 2016 at 9:51 AM, Jackson Cavett <jcavett_at_mail.bradley.edu>
wrote:

> I gave it a try and ran into the same problem. I didn't do any constraint
> scaling, and the following are the harmonic restraint parameters I used:
>
> constantforce yes
> consforcefile alphaV_equil_harm_restraint.ref
>
> constraints on
> consref alphaV_equil.pdb
> conskfile alphaV_equil_harm_restraint.ref
> conskcol B
> constraintScaling 1.0
>
>
> Please let me know if I have anything setup wrong or if there is a change
> that might help. I have attached a screenshot of what the output looks like
> after about 4 picoseconds. For this simulation I used a force of 1
> kcal/mol-angstrom^2, which I don't think is huge. Hopefully that picture
> can help make sense of what is going wrong.
>
> Thanks Again,
>
> Jack
>
>
>
>
>
> On Fri, Aug 12, 2016 at 3:56 PM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov>
> wrote:
>
>> Ok. Have you tried using constraints rather than fixing the atom
>> positions? (http://www.ks.uiuc.edu/Research/namd/2.11/ug/node27.html#SE
>> CTION00086200000000000000) In my limited experience with SMD, fixing
>> atoms causes problems when particles get close to the fixed atoms, and you
>> get failures of those specific atoms that venture near the brick wall. Your
>> minimization is probably sufficient. As a general rule, 1000 steps is
>> normally enough to start a unequilibrated simulation, and 2000 might be
>> needed for tough cases. Anything more than that generally doesn't help me,
>> and instead instabilities are the result of me doing something dumb (like
>> piercing a ring or putting two atoms exactly on top of one another). But,
>> since this came from an existing trajectory, its probably a fixed atom
>> problem.
>>
>> -Josh
>>
>> On 08/12/2016 02:21 PM, Jackson Cavett wrote:
>> Hi Josh,
>>
>> Other research using the same molecule have used forces around these
>> numbers, so that is where I got started. In my reference file, I set the
>> protein residues 1-440 as the SMD atoms and set protein residues 1500-1700
>> fixed. The SMD direction was the vector from the fixed atoms to the SMD
>> atoms. In the configuration file, I changed the dcd frequency and other
>> frequencies to 1 and ran a very short simulation. My real simulations will
>> hopefully use a frequency of 500 and run for somewhere between 1 and 5
>> nanoseconds. My configuration file looks as follows:
>>
>> #############################################################
>> ## JOB DESCRIPTION ##
>> #############################################################
>>
>> # Constant Force Pulling for fully equilibrated alphaVbeta3
>> # 0.5 nanoseconds with a force of 10 (about 700 pN)
>>
>>
>> #############################################################
>> ## ADJUSTABLE PARAMETERS ##
>> #############################################################
>>
>> structure alphaVbeta3_wb.psf
>> coordinates alphaV_full_equil.pdb
>> outputName first_fully_equilibrated_alphaV
>>
>> set temperature 310
>>
>> # Continuing a job from the restart files
>> #if {0} {
>> #set inputname myinput
>> #binCoordinates $inputname.restart.coor
>> #binVelocities $inputname.restart.vel ;# remove the "temperature"
>> entry if you use this!
>> #extendedSystem $inputname.xsc
>> #}
>>
>> firsttimestep 0
>>
>>
>> #############################################################
>> ## SIMULATION PARAMETERS ##
>> #############################################################
>>
>> # Input
>> paraTypeCharmm on
>> parameters ../par_all27_prot_lipid.inp
>>
>> # NOTE: Do not set the initial velocity temperature if you
>> # have also specified a .vel restart file!
>> temperature $temperature
>>
>>
>> # Periodic Boundary conditions
>> # NOTE: Do not set the periodic cell basis if you have also
>> # specified an .xsc restart file!
>> if {1} {
>> cellBasisVector1 200.0 0 0
>> cellBasisVector2 0 200.0 0
>> cellBasisVector3 0 0 250.0
>> cellOrigin -16.3 33.85 38.46
>> }
>> wrapWater on
>> wrapAll on
>>
>>
>> # Force-Field Parameters
>> exclude scaled1-4
>> 1-4scaling 1.0
>> cutoff 12.0
>> switching on
>> switchdist 10.0
>> pairlistdist 14.0
>>
>>
>> # Integrator Parameters
>> timestep 2.0 ;# 2fs/step
>> rigidBonds all ;# needed for 2fs steps
>> nonbondedFreq 1
>> fullElectFrequency 2
>> stepspercycle 10
>>
>>
>> #PME (for full-system periodic electrostatics)
>> if {0} {
>> PME yes
>> PMEGridSpacing 1.0
>>
>> #manual grid definition
>> #PMEGridSizeX 32
>> #PMEGridSizeY 32
>> #PMEGridSizeZ 64
>> #}
>>
>>
>> # Constant Temperature Control
>> langevin off ;# do langevin dynamics
>> langevinDamping 1 ;# damping coefficient (gamma) of 5/ps
>> langevinTemp $temperature
>> langevinHydrogen no ;# don't couple langevin bath to hydrogens
>>
>>
>> # Constant Pressure Control (variable volume)
>> if {0} {
>> useGroupPressure yes ;# needed for 2fs steps
>> useFlexibleCell no ;# no for water box, yes for membrane
>> useConstantArea no ;# no for water box, yes for membrane
>>
>> langevinPiston on
>> langevinPistonTarget 1.01325 ;# in bar -> 1 atm
>> langevinPistonPeriod 100.0
>> langevinPistonDecay 50.0
>> langevinPistonTemp $temperature
>> }
>>
>>
>> restartfreq 1 ;# 500steps = every 1ps
>> dcdfreq 1
>> xstFreq 1
>> outputEnergies 1
>> outputPressure 100
>>
>>
>> # Fixed Atoms Constraint (set PDB beta-column to 1)
>> if {1} {
>> fixedAtoms on
>> fixedAtomsFile first_sim_attempt.ref
>> fixedAtomsCol B
>> }
>>
>>
>> # IMD Settings (can view sim in VMD)
>> if {0} {
>> IMDon on
>> IMDport 3000 ;# port number (enter it in VMD)
>> IMDfreq 1 ;# send every 1 frame
>> IMDwait no ;# wait for VMD to connect before running?
>> }
>>
>>
>> #############################################################
>> ## EXTRA PARAMETERS ##
>> #############################################################
>>
>> # Put here any custom parameters that are specific to
>> # this job (e.g., SMD, TclForces, etc...)
>>
>> constantforce yes
>> consforcefile first_sim_attempt.ref
>>
>>
>> #############################################################
>> ## EXECUTION SCRIPT ##
>> #############################################################
>>
>> #Minimization
>> if {0} {
>> minimize 10000
>> reinitvels $temperature
>> }
>>
>> run 250 ;# 0.5ns
>>
>> On Fri, Aug 12, 2016 at 2:45 PM, Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov
>> <mailto:Joshua.Vermaas_at_nrel.gov>> wrote:
>> Hi Jack,
>>
>> Usually, the SMD forces should be high enough to force movement, but not
>> high enough that everything else stays still by comparison, especially over
>> several timesteps. Odds are you are driving your system far too hard, and
>> crashing as a result. What does your NAMD configuration file look like?
>>
>> -Josh
>>
>> On 08/12/2016 01:31 PM, Jackson Cavett wrote:
>> Hello All,
>>
>> I have been trying to equilibrate and minimize integrin alphaVbeta3 in a
>> waterbox. I finished up a long equilibration simulation (about 40 ns), but
>> am running into the "ERROR: Constraint failure in RATTLE algorithm for atom
>> 994!" when I try to run a constant-force simulation. After reading up on
>> this error, I changed the dcd frequency to every timestep and looked at the
>> dcd file in VMD. When looking at the dcd file, it looks as though the smd
>> atoms being pulled (by about 700 pN) move while every other atom stays
>> fixed. It seems to me (though I am far from an expert) that the energy has
>> not reached a constant value and therefore I would need more minimization
>> steps. Before the equilibration, I minimized the system with three
>> different minimzation simulations which were 10,000 timesteps each. My
>> system has about 170,000 atoms, which is very large. Would anyone know how
>> many minimization steps this system would need? Also, I could be on the
>> wrong track with what my problem is, so if anyone has any other ideas that
>> would be great. Thank you all for your help!
>>
>> Thanks,
>>
>> Jack
>>
>>
>>
>>
>