Re: Trajectory analysis

From: Philip Blood (philb_at_hec.utah.edu)
Date: Tue Oct 17 2006 - 18:42:56 CDT

Hi,

Was there a resolution to this problem? I have run across a similar
problem doing trajectory analysis in NAMD2.6b1 (currently running on
NCSA IA-64). I found that when analyzing an NPT trajectory using the Tcl
functions and reading in a dcd file, it appears that NAMD does not
update the cell dimensions when it calculates the forces. This is
manifested in a number of ways: first, I expected that when analyzing
the dcd I would not need to worry about entering the precise cell
dimensions into the NAMD configuration file since it would read them
from the dcd. This, however, resulted in immediate constraint
failures. So I entered the exact cell dimensions in the configuration
file (obtained from the .xst file from the simulation that generated the
dcd), and the simulation ran for a few steps and then failed with
constraint failures. In the initial step the energy looks perfect. If
I turn off rigidBonds, I see the same problem mentioned below: the
analysis runs, but eventually the energies become huge.

I used the getcell command in a tclBC script to check the cell
dimensions, and NAMD is reading them correctly, but for some reason is
not using them when it calculates the force--I believe it just uses the
original ones that I input directly into the configuration file for the
purpose of calculating the forces and energy. In addition, my colleague
has done something similar and obtained the same result, however, if she
then ran the same analysis on an NVT simulation trajectory everything
was perfect.

Here is that script and the output that shows it reading the right cell
dimensions but then terminating with constraint failures. Again, if I
remove the rigidBonds it continues to run, but the energies blow up.

Thanks,
Phil Blood

source source.file

#firsttimestep $first
#numsteps $numsteps

temperature 310K
seed 40277

coordinates $dir0/$name.pdb
#bincoordinates $dir0/$name.$id0.restart.coor
#binvelocities $dir0/$name.$id0.restart.vel

structure $dir0/$name.xplor.psf

parameters $dir0/par_dops.inp
paraTypeCharmm on

cellBasisVector1 459.122 0.0 0.0
cellBasisVector2 0.0 93.9624 0.0
cellBasisVector3 0.0 0.0 167.397
#extendedSystem $dir0/$name.$id0.restart.xsc
XSTfile $scrdir/$name.$id.xst
XSTfreq 500
wrapAll on

exclude scaled1-4
1-4scaling 1.0
switching on
switchdist 10.0
cutoff 12.0
pairlistdist 13.5
stepspercycle 20
#stepspercycle 1
rigidBonds all
rigidTolerance 1.0e-8

#PME options
PME yes
PMETolerance 1.0e-6
PMEInterpOrder 6
PMEGridSizeX 486
PMEGridSizeY 100
PMEGridSizeZ 150

# pressure control
useGroupPressure yes
useFlexibleCell yes
useConstantRatio no
useConstantArea no
LangevinPiston on
LangevinPistonTarget 1.01325
LangevinPistonPeriod 2000
LangevinPistonDecay 2000
LangevinPistonTemp 310

# temperature control
langevin on
langevinTemp 310
langevinDamping 0.5

fullElectFrequency 1
timestep 2.0

# output params
outputname $scrdir/$name.$id
#outputEnergies 500
#outputMomenta 500
#outputPressure 500

# outputTiming
#binaryoutput no
#DCDfile $dcddir/$name.$id.dcd
#DCDfreq 500
DCDUnitCell yes
#velDCDfile $dcddir/$name.$id.veldcd
#velDCDfreq 500
#restartname $scrdir/$name.$id.restart
#restartfreq 500

tclBC on

#files for input and output
set inputdcd $dir0/$name.$id.dcd
set ftmp $scrdir/$name.fout.$id.dat
set f_outfile [open $ftmp w]

tclBCScript {
 proc calcforces {step unique} {

  set mycell [getcell]
  print $mycell

  while {[nextatom]} {

  }
 }
}

set ts 500
firstTimestep $ts
set maxstep 50000
coorfile open dcd $inputdcd

while { ![coorfile read] && $ts < $maxstep } {
# Compute energies and forces, but don't try to move the atoms.
  run 0
  puts "reading timestep $ts"
  incr ts 500
  firstTimestep $ts
}

run 0

close $f_outfile

Info: NAMD 2.6b1 for Linux-ia64-MPI-GM
Info:
Info: Please visit http://www.ks.uiuc.edu/Research/namd/
Info: and send feedback or bug reports to namd_at_ks.uiuc.edu
Info:
Info: Please cite Kale et al., J. Comp. Phys. 151:283-312 (1999)
Info: in all publications reporting results obtained with NAMD.
Info:
Info: Based on Charm++/Converse 50900 for mpi-linux-ia64
Info: Built Fri Jul 29 11:37:35 CDT 2005 by jphillip on tg-login2
Info: Sending usage information to NAMD developers via UDP. Sent data is:
Info: 1 NAMD 2.6b1 Linux-ia64-MPI-GM 4 tg-c062 philipblood
Info: Running on 4 processors.
Info: 20819 kB of memory in use.
Info: Changed directory to
/home/philipblood/work/bar/force-match/nbar_dops.10x2.fm
Info: Configuration file is d.namd
TCL: Suspending until startup complete.
Info: SIMULATION PARAMETERS:
Info: TIMESTEP 2
Info: NUMBER OF STEPS 500
Info: STEPS PER CYCLE 20
Info: PERIODIC CELL BASIS 1 459.122 0 0
Info: PERIODIC CELL BASIS 2 0 93.9624 0
Info: PERIODIC CELL BASIS 3 0 0 167.397
Info: PERIODIC CELL CENTER 0 0 0
Info: WRAPPING ALL CLUSTERS AROUND PERIODIC BOUNDARIES ON OUTPUT.
Info: LOAD BALANCE STRATEGY Other
Info: LDB PERIOD 4000 steps
Info: FIRST LDB TIMESTEP 100
Info: LDB BACKGROUND SCALING 1
Info: HOM BACKGROUND SCALING 1
Info: PME BACKGROUND SCALING 1
Info: MAX SELF PARTITIONS 50
Info: MAX PAIR PARTITIONS 20
Info: SELF PARTITION ATOMS 125
Info: PAIR PARTITION ATOMS 200
Info: PAIR2 PARTITION ATOMS 400
Info: MIN ATOMS PER PATCH 100
Info: INITIAL TEMPERATURE 310
Info: CENTER OF MASS MOVING? NO
Info: DIELECTRIC 1
Info: EXCLUDE SCALED ONE-FOUR
Info: 1-4 SCALE FACTOR 1
Info: NO DCD TRAJECTORY OUTPUT
Info: XST FILENAME
/disks/scratchgpfs1/philipblood/nbar_dops.10x2.fm/339492.d179/nbar_dops.10x2.fm.d179.xst
Info: XST FREQUENCY 500
Info: NO VELOCITY DCD OUTPUT
Info: OUTPUT FILENAME
/disks/scratchgpfs1/philipblood/nbar_dops.10x2.fm/339492.d179/nbar_dops.10x2.fm.d179
Info: BINARY OUTPUT FILES WILL BE USED
Info: NO RESTART FILE
Info: SWITCHING ACTIVE
Info: SWITCHING ON 10
Info: SWITCHING OFF 12
Info: PAIRLIST DISTANCE 13.5
Info: PAIRLIST SHRINK RATE 0.01
Info: PAIRLIST GROW RATE 0.01
Info: PAIRLIST TRIGGER 0.3
Info: PAIRLISTS PER CYCLE 2
Info: PAIRLISTS ENABLED
Info: MARGIN 0.96
Info: HYDROGEN GROUP CUTOFF 2.5
Info: PATCH DIMENSION 16.96
Info: TCL BOUNDARY FORCES ACTIVE
Info: TCL BOUNDARY FORCES SCRIPT
proc calcforces {step unique} {

set mycell [getcell]
print $mycell

while {[nextatom]} {

}

}

Info: TCL BOUNDARY FORCES ARGS
Info: LANGEVIN DYNAMICS ACTIVE
Info: LANGEVIN TEMPERATURE 310
Info: LANGEVIN DAMPING COEFFICIENT IS 0.5 INVERSE PS
Info: LANGEVIN DYNAMICS APPLIED TO HYDROGENS
Info: LANGEVIN PISTON PRESSURE CONTROL ACTIVE
Info: TARGET PRESSURE IS 1.01325 BAR
Info: OSCILLATION PERIOD IS 2000 FS
Info: DECAY TIME IS 2000 FS
Info: PISTON TEMPERATURE IS 310 K
Info: PRESSURE CONTROL IS GROUP-BASED
Info: INITIAL STRAIN RATE IS 0 0 0
Info: CELL FLUCTUATION IS ANISOTROPIC
Info: PARTICLE MESH EWALD (PME) ACTIVE
Info: PME TOLERANCE 1e-06
Info: PME EWALD COEFFICIENT 0.257952
Info: PME INTERPOLATION ORDER 6
Info: PME GRID DIMENSIONS 486 100 150
Info: Attempting to read FFTW data from
FFTW_NAMD_2.6b1_Linux-ia64-MPI-GM.txt
Info: Optimizing 6 FFT steps. 1... 2... 3... 4... 5... 6... Done.
Info: Writing FFTW data to FFTW_NAMD_2.6b1_Linux-ia64-MPI-GM.txt
Info: FULL ELECTROSTATIC EVALUATION FREQUENCY 1
Info: USING VERLET I (r-RESPA) MTS SCHEME.
Info: C1 SPLITTING OF LONG RANGE ELECTROSTATICS
Info: PLACING ATOMS IN PATCHES BY HYDROGEN GROUPS
Info: RIGID BONDS TO HYDROGEN : ALL
Info: ERROR TOLERANCE : 1e-08
Info: MAX ITERATIONS : 100
Info: RIGID WATER USING SETTLE ALGORITHM
Info: RANDOM NUMBER SEED 40277
Info: USE HYDROGEN BONDS? NO
Info: COORDINATE PDB
/home/philipblood/work/bar/force-match/nbar_dops.10x2.fm/nbar_dops.10x2.fm.pdb
Info: STRUCTURE FILE
/home/philipblood/work/bar/force-match/nbar_dops.10x2.fm/nbar_dops.10x2.fm.xplor.psf
Info: PARAMETER file: CHARMM format!
Info: PARAMETERS
/home/philipblood/work/bar/force-match/nbar_dops.10x2.fm/par_dops.inp
Info: USING ARITHMETIC MEAN TO COMBINE L-J SIGMA PARAMETERS
Info: FIRST TIMESTEP 500
Info: SUMMARY OF PARAMETERS:
Info: 181 BONDS
Info: 454 ANGLES
Info: 570 DIHEDRAL
Info: 46 IMPROPER
Info: 114 VDW
Info: 0 VDW_PAIRS
Warning: Ignored 177536 bonds with zero force constants.
Warning: Will get H-H distance in rigid H2O from H-O-H angle.
Info: ****************************
Info: STRUCTURE SUMMARY:
Info: 737674 ATOMS
Info: 557330 BONDS
Info: 562792 ANGLES
Info: 534758 DIHEDRALS
Info: 4718 IMPROPERS
Info: 0 EXCLUSIONS
Info: 654538 RIGID BONDS
Info: 1558484 DEGREES OF FREEDOM
Info: 260672 HYDROGEN GROUPS
Info: TOTAL MASS = 4.42599e+06 amu
Info: TOTAL CHARGE = 0.000193976 e
Info: *****************************
Info: Entering startup phase 0 with 208622 kB of memory in use.
Info: Entering startup phase 1 with 208622 kB of memory in use.
Info: Entering startup phase 2 with 327085 kB of memory in use.
Info: Entering startup phase 3 with 332861 kB of memory in use.
Info: PATCH GRID IS 27 (PERIODIC) BY 5 (PERIODIC) BY 9 (PERIODIC)
Info: REMOVING COM VELOCITY -0.00906089 0.00136749 -0.00782743
Info: LARGEST PATCH (674) HAS 1093 ATOMS
Info: Entering startup phase 4 with 428219 kB of memory in use.
Info: PME using 4 and 4 processors for FFT and reciprocal sum.
Info: PME GRID LOCATIONS: 0 1 2 3
Info: PME TRANS LOCATIONS: 0 1 2 3
Info: Optimizing 4 FFT steps. 1... 2... 3... 4... Done.
Info: Entering startup phase 5 with 442739 kB of memory in use.
Info: Entering startup phase 6 with 376170 kB of memory in use.
Measuring processor speeds... Done.
Info: Entering startup phase 7 with 376263 kB of memory in use.
Info: NONBONDED TABLE R-SQUARED SPACING: 0.0625
Info: NONBONDED TABLE SIZE: 769 POINTS
Info: ABSOLUTE IMPRECISION IN FAST TABLE FORCE: 2.64698e-23 AT 11.9138
Info: RELATIVE IMPRECISION IN FAST TABLE FORCE: 1.40436e-16 AT 11.9138
Info: ABSOLUTE IMPRECISION IN SCOR TABLE FORCE: 4.10282e-22 AT 11.9138
Info: RELATIVE IMPRECISION IN SCOR TABLE FORCE: 3.56912e-15 AT 11.9138
Info: ABSOLUTE IMPRECISION IN VDWA TABLE ENERGY: 1.75 AT 0.0441942
Info: RELATIVE IMPRECISION IN VDWA TABLE ENERGY: 6.75473e-14 AT 11.9138
Info: ABSOLUTE IMPRECISION IN VDWA TABLE FORCE: 19968 AT 0.0441942
Info: RELATIVE IMPRECISION IN VDWA TABLE FORCE: 6.16499e-15 AT 0.0441942
Info: ABSOLUTE IMPRECISION IN VDWB TABLE ENERGY: 4.1359e-25 AT 11.8295
Info: RELATIVE IMPRECISION IN VDWB TABLE ENERGY: 7.56853e-15 AT 11.9138
Info: ABSOLUTE IMPRECISION IN VDWB TABLE FORCE: 3.87741e-26 AT 11.9138
Info: RELATIVE IMPRECISION IN VDWB TABLE FORCE: 5.97409e-16 AT 11.9138
Info: Entering startup phase 8 with 409365 kB of memory in use.
Info: Finished startup with 450695 kB of memory in use.
Info: Coordinate file
/home/philipblood/work/bar/force-match/nbar_dops.10x2.fm/nbar_dops.10x2.fm.d179.dcd
opened for reading.
Info: Reading timestep from file.
Info: Updating unit cell from timestep.
Info: REMOVING COM VELOCITY -0.00906089 0.00136749 -0.00782743
TCL: Running for 0 steps
TCL: {0.0 0.0 0.0} {459.122253418 0.0 0.0} {-0.0 93.9624481201 0.0}
{-0.0 -0.0 167.396759033}
TCL: {0.0 0.0 0.0} {459.122253418 0.0 0.0} {-0.0 93.9624481201 0.0}
{-0.0 -0.0 167.396759033}
TCL: {0.0 0.0 0.0} {459.122253418 0.0 0.0} {-0.0 93.9624481201 0.0}
{-0.0 -0.0 167.396759033}
TCL: {0.0 0.0 0.0} {459.122253418 0.0 0.0} {-0.0 93.9624481201 0.0}
{-0.0 -0.0 167.396759033}
ETITLE: TS BOND ANGLE DIHED
IMPRP ELECT VDW BOUNDARY
MISC KINETIC TOTAL TEMP
TOTAL2 TOTAL3 TEMPAVG PRESSURE
GPRESSURE VOLUME PRESSAVG GPRESSAVG

ENERGY: 500 24056.0398 122674.8128 54536.6888
1200.1142 -2235410.8602 113094.2910 0.0000
0.0000 480348.6491 -1439500.2646 310.2020 -1436969.0327
-1436969.0327 310.2020 -26.5192 -25.0174
7221538.1875 -26.5192 -25.0174

OPENING EXTENDED SYSTEM TRAJECTORY FILE
reading timestep 500
TCL: Setting parameter firstTimestep to 1000
Info: Reading timestep from file.
Info: Updating unit cell from timestep.
Info: REMOVING COM VELOCITY -0.00906089 0.00136749 -0.00782743
TCL: Running for 0 steps
TCL: {0.0 0.0 0.0} {458.889312744 0.0 0.0} {-0.0 94.0012283325 0.0}
{-0.0 -0.0 167.454818726}
TCL: {0.0 0.0 0.0} {458.889312744 0.0 0.0} {-0.0 94.0012283325 0.0}
{-0.0 -0.0 167.454818726}
TCL: {0.0 0.0 0.0} {458.889312744 0.0 0.0} {-0.0 94.0012283325 0.0}
{-0.0 -0.0 167.454818726}
TCL: {0.0 0.0 0.0} {458.889312744 0.0 0.0} {-0.0 94.0012283325 0.0}
{-0.0 -0.0 167.454818726}
ETITLE: TS BOND ANGLE DIHED
IMPRP ELECT VDW BOUNDARY
MISC KINETIC TOTAL TEMP
TOTAL2 TOTAL3 TEMPAVG PRESSURE
GPRESSURE VOLUME PRESSAVG GPRESSAVG

ENERGY: 1000 24673.4115 123142.5547 54331.6447
1211.3140 -2237759.5570 115041.8421 0.0000
0.0000 480753.6046 -1438605.1854 310.4635 -1436057.1778
-1436120.0866 310.4635 16.4442 17.2117
7223357.6970 16.4442 17.2117

reading timestep 1000
TCL: Setting parameter firstTimestep to 1500
Info: Reading timestep from file.
Info: Updating unit cell from timestep.
Info: REMOVING COM VELOCITY -0.00906089 0.00136749 -0.00782743
TCL: Running for 0 steps
TCL: {0.0 0.0 0.0} {458.79006958 0.0 0.0} {-0.0 93.9584503174 0.0} {-0.0
-0.0 167.486663818}
TCL: {0.0 0.0 0.0} {458.79006958 0.0 0.0} {-0.0 93.9584503174 0.0} {-0.0
-0.0 167.486663818}
TCL: {0.0 0.0 0.0} {458.79006958 0.0 0.0} {-0.0 93.9584503174 0.0} {-0.0
-0.0 167.486663818}
TCL: {0.0 0.0 0.0} {458.79006958 0.0 0.0} {-0.0 93.9584503174 0.0} {-0.0
-0.0 167.486663818}
ETITLE: TS BOND ANGLE DIHED
IMPRP ELECT VDW BOUNDARY
MISC KINETIC TOTAL TEMP
TOTAL2 TOTAL3 TEMPAVG PRESSURE
GPRESSURE VOLUME PRESSAVG GPRESSAVG

ENERGY: 1500 25155.0185 123729.6189 54457.3885
1228.1731 -2237642.8195 114782.3239 0.0000
0.0000 480551.7452 -1437738.5514 310.3332 -1435167.5579
-1435312.7320 310.3332 159.5641 159.0790
7219881.7776 159.5641 159.0790

reading timestep 1500
TCL: Setting parameter firstTimestep to 2000
Info: Reading timestep from file.
Info: Updating unit cell from timestep.
Info: REMOVING COM VELOCITY -0.00906089 0.00136749 -0.00782743
TCL: Running for 0 steps
TCL: {0.0 0.0 0.0} {459.304840088 0.0 0.0} {-0.0 93.8790893555 0.0}
{-0.0 -0.0 167.392944336}
TCL: {0.0 0.0 0.0} {459.304840088 0.0 0.0} {-0.0 93.8790893555 0.0}
{-0.0 -0.0 167.392944336}
TCL: {0.0 0.0 0.0} {459.304840088 0.0 0.0} {-0.0 93.8790893555 0.0}
{-0.0 -0.0 167.392944336}
TCL: {0.0 0.0 0.0} {459.304840088 0.0 0.0} {-0.0 93.8790893555 0.0}
{-0.0 -0.0 167.392944336}
ETITLE: TS BOND ANGLE DIHED
IMPRP ELECT VDW BOUNDARY
MISC KINETIC TOTAL TEMP
TOTAL2 TOTAL3 TEMPAVG PRESSURE
GPRESSURE VOLUME PRESSAVG GPRESSAVG

ENERGY: 2000 24868.6147 123753.2825 54229.5853
1209.7959 -2235748.9915 113579.3721 0.0000
0.0000 480274.4063 -1437833.9347 310.1541 -1435275.6421
-1435364.1144 310.1541 34.7683 32.3176
7217836.4747 34.7683 32.3176

reading timestep 2000
TCL: Setting parameter firstTimestep to 2500
Info: Reading timestep from file.
Info: Updating unit cell from timestep.
Info: REMOVING COM VELOCITY -0.00906089 0.00136749 -0.00782743
TCL: Running for 0 steps
TCL: {0.0 0.0 0.0} {459.766082764 0.0 0.0} {-0.0 93.8565673828 0.0}
{-0.0 -0.0 167.408905029}
TCL: {0.0 0.0 0.0} {459.766082764 0.0 0.0} {-0.0 93.8565673828 0.0}
{-0.0 -0.0 167.408905029}
TCL: {0.0 0.0 0.0} {459.766082764 0.0 0.0} {-0.0 93.8565673828 0.0}
{-0.0 -0.0 167.408905029}
TCL: {0.0 0.0 0.0} {459.766082764 0.0 0.0} {-0.0 93.8565673828 0.0}
{-0.0 -0.0 167.408905029}
ERROR: Constraint failure in RATTLE algorithm for atom 160079!
ERROR: Constraint failure; simulation has become unstable.
ERROR: Constraint failure in RATTLE algorithm for atom 339349!
ERROR: Constraint failure; simulation has become unstable.
ERROR: Constraint failure in RATTLE algorithm for atom 336289!
ERROR: Constraint failure; simulation has become unstable.
ERROR: Exiting prematurely.
==========================================
WallClock: 113.538277 CPUTime: 113.538277 Memory: 662465 kB
End of program

Peter Freddolino wrote:
> Hi Luca,
> could you let me know what versions of VMD and NAMD you're using, and
> what OS? And can I verify your procedure: You loaded your psf and dcd
> into vmd and then ran namdenergy? How many frames did your dcd have? It
> looks like there's only one frame of namdenergy output...
>
> Peter
>
> Luca Bellucci wrote:
>
>> Hi all,
>> I have dcd file named file.dcd for a MD of this kind:
>> NPT PME with these option
>> dcdfreq 150
>> xstFreq 150
>> wrapAll on
>> dcdUnitCell yes
>> Then I have this script (arranged by NAMDEnergy VMD script) for
>> post-analysis:
>>
>> \ structure file.psf
>> \ paraTypeCharmm on
>> \ parameters par_all27_prot_lipid_na.inp
>> \ numsteps 1
>> \ exclude scaled1-4
>> \ outputname namd-temp
>> \ temperature 0
>> \ cutoff 12
>> \ switchdist 10
>> \ pairInteraction on
>> \ pairInteractionGroup1 1
>> \ pairInteractionFile namd-temp.pdb
>> \ pairInteractionGroup2 2
>> \ coordinates namd-temp.pdb
>> \ run 0
>> \ set ts 0
>> \ coorfile open dcd namd-temp.dcd
>> \ while { ![coorfile read] } {
>> \ firstTimestep 0
>> \ run 0
>> \ incr ts 1
>> \ }
>> \ coorfile close
>>
>> But the output for energy is wrong. For example:
>> ELECT VDW
>> 1427.4544 99999999.9999
>> File namd-temp.dcd generated by
>> NAMDEnergy script from my original file.dcd. Is this a problem?
>> Have you any solution?
>>
>> Thanks
>> Luca
>>

-- 
Philip D. Blood
Research Assistant
Center for Biophysical Modeling and Simulation
Dept. of Bioengineering, University of Utah
315 S. 1400 E., RM. 2020
Salt Lake City, UT 84112-0850
Phone: (801) 581-8606; Fax: (801) 581-4353

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:44:04 CST