I am simulating a protein-protein complex and have placed the complex in a minimal water box with 15 Angstroms of distance from the end of the box. Now that I have completed the simulation, I ran RMSD analysis using the following script:

-------------------------------------------------------------------------------------

# wrap a trajectory to avoid RMSD calculation errors

pbc wrap -centersel "protein" -center com -compound residue -all

## Align to first frame

# the frame being compared

set reference [atomselect top "protein and backbone" frame 0]

# the frame to be compared

set compare [atomselect top "protein and backbone"]

#get the number of frames of the trajectory

set num_steps [molinfo top get numframes]

# get the correct num as the trajectory start from zero

set outfile [open rmsd.dat w]

for {set frame 0} {$frame < $num_steps} {incr frame} {

# get the correct frame

$compare frame $frame

# compute the 4*4 matrix transformation that takes one set of coordinates onto the other

set trans_mat [measure fit $compare $reference]

# do the alignment

$compare move $trans_mat

# compute the RMSD

set rmsd [measure rmsd $compare $reference ]

# print the RMSD

puts $outfile "$frame $rmsd"

}

close $outfile

------------------------------------------------------------------------------------

When checking the results, the RMSD starts off at a reasonable 1-2 Angstroms, but then rapidly and randomly increases to over 10 Angstroms (image attached).

Watching the simulation on VMD, I can see that initially the protein is located in the center of the water box. However, as the simulation proceeds, the protein goes to the edge of the box and then appears on the other side. This seems to coincide with the very strange RMSD results.

Is there a way to ensure that the protein does not leave the center of the water box? Or do I need to modify my RMSD script in some way to prevent these strange results?

Here is the MD.conf file used for the simulation:

-----------------------------------------------------------------------------------

cutoff 12.0

pairlistdist 14.0

switching on

switchdist 10.0

PME on

PMEGridspacing 1

wrapAll on

wrapWater on

############################################################################

#cr

#cr (C) Copyright 1995-2009 The Board of Trustees of the

#cr University of Illinois

#cr All Rights Reserved

#cr

############################################################################

############################################################################

# RCS INFORMATION:

#

# $RCSfile: MD.conf,v $

# $Author: jribeiro $ $Locker: $ $State: Exp $

# $Revision: 1.2 $ $Date: 2017/05/10 19:03:08 $

#

############################################################################

##START HERE###

##Simulation Template##

# Simulation conditions

coordinates C3d_G321D_100ns_Test1_QwikMD.pdb

structure C3d_G321D_100ns_Test1_QwikMD.psf

binCoordinates Equilibration.4.restart.coor

binVelocities Equilibration.4.restart.vel

extendedSystem Equilibration.4.restart.xsc

# Simulation conditions

#temperature 0

# Harmonic constraints

constraints off

consref qwikmdTemp_constraints.pdb

conskfile qwikmdTemp_constraints.pdb

constraintScaling 2

consexp 2

conskcol B

# Output Parameters

binaryoutput no

outputname MD

outputenergies 40

outputtiming 40

outputpressure 40

binaryrestart yes

dcdfile MD.dcd

dcdfreq 1000

XSTFreq 1000

restartfreq 1000

restartname MD.restart

# Thermostat Parameters

langevin on

langevintemp 310

langevinHydrogen off

langevindamping 1

# Barostat Parameters

usegrouppressure yes

useflexiblecell no

useConstantArea no

langevinpiston on

langevinpistontarget 1.01325

langevinpistonperiod 200

langevinpistondecay 100

langevinpistontemp 310

# Integrator Parameters

timestep 2

firstTimestep 0

fullElectFrequency 2

nonbondedfreq 1

stepspercycle 10

# Force Field Parameters

paratypecharmm on

parameters toppar_water_ions_namd.str

parameters toppar_all36_carb_glycopeptide.str

parameters par_all36_lipid.prm

parameters par_all36_na.prm

parameters par_all36_prot.prm

parameters par_all36_carb.prm

parameters par_all36_cgenff.prm

exclude scaled1-4

1-4scaling 1.0

rigidbonds all

#Implicit Solvent Parameters

gbis off

alphaCutoff 14.0

ionConcentration 0.15

# Script

run 50000000

set file [open MD.check w+]

set done 1

if {[file exists MD.restart.coor] != 1 || [file exists MD.restart.vel] != 1 || [file exists MD.restart.xsc] != 1 } {

set done 0

}

if {$done == 1} {

puts $file "DONE"

flush $file

close $file

} else {

puts $file "One or more files failed to be written"

flush $file

close $file

}

--------------------------------------------------------------------------------------