VMD-L Mailing List

From: Vermaas, Joshua (Joshua.Vermaas_at_nrel.gov)
Date: Sat Aug 13 2016 - 15:37:45 CDT

Hi Jack,

That is one wickedly extended protein. Was it always like that? It looks like the piece in the lower right blew out of the waterbox at quite a high velocity, so it isn't surprising to me that this simulation crashed.

-Josh

On 08/13/2016 08:51 AM, Jackson Cavett 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<mailto: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#SECTION00086200000000000000) 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><mailto: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 t
hat would be great. Thank you all for your help!

Thanks,

Jack