Runaway positive feedback when minimizing large systems (was: strange glitches when performing energy minimization in MDFF simulations)

From: Tristan Croll (tristan.croll_at_qut.edu.au)
Date: Tue Feb 11 2014 - 20:07:59 CST

May I ask anyone who's interested to perform the following experiment, please?

Take a large, equilibrated system of your choice (e.g. the ATPase benchmark). Forcefield doesn't seem to matter. Choose an asparagine residue somewhere near the center, and move OD1 right up against CG. Then set up a minimization AutoIMD simulation as "same residue as within 10 of [your chosen residue]". If it behaves the same way as mine, this will very quickly snap back to a normal arrangement. Discard that simulation, and change the text to "same residue as within 50 of...". In this case, you should see behaviour reminiscent of the videos below.

What I think is happening is that as the dimensions of the simulation box increase, the coordinate grid becomes too coarse for the minimization function to work properly. Past a certain scale (the dimensions of the system in which I first noticed a problem are {153.914001 175.417007 186.804001}, even structures that are perhaps statistical outliers but perfectly normal results of equilibrium simulations can be enough to throw it into a fit.

Cheers,

Tristan

From: owner-namd-l_at_ks.uiuc.edu [mailto:owner-namd-l_at_ks.uiuc.edu] On Behalf Of Tristan Croll
Sent: Wednesday, 12 February 2014 7:33 AM
To: namd-l_at_ks.uiuc.edu
Subject: Re: namd-l: RE: Strange glitches when performing energy minimization in MDFF simulations

I've tracked the problem down some more. Firstly, it's not hardware (or software version) related - I've tried it on the cluster, and on my own workstation with 64-bit NAMD 2.8 and 2.9, and with NAMD 2.9-CUDA. Secondly, while the large spikes in VdW energy occur later in the minimisation, the problem starts at the first time step. Thirdly, it's reproducible - for a particular starting configuration the problem will always arise in the same place - one (perfectly normal looking, thoroughly equilibrated) residue will go haywire. The energy spikes I noticed seem to arise when this leads to two atoms attempting to occupy the same point in space. Fourthly, it has nothing to do with MDFF - the behaviour arises whether or not grid forces are on. What's particularly confusing is that with the same PSF file but different coordinates (I.e. taken from different points in the trajectory) the behaviour will either disappear entirely or show up in an entirely different location.

In most cases, it's all over by step 1000, when I would normally have NAMD write the first frame of the trajectory. In order to get a better picture of what's going on, I set the dcd output instead to 10 time step intervals. These movies show what I saw:

https://www.dropbox.com/s/oy8xgcbgy1e2k8g/restart_1.mpg

https://www.dropbox.com/s/rqzxmta311vxgs9/restart_2.mpg

The construct is essentially the same in each of these cases. The only difference is that in the one showing the histidine residue the construct has been rebuilt after reassigning HIS protonation states. This one is a straightforward restart after >10ns of uneventful equilibration; the other is starting from a rewritten PDB file after performing an interactive simulation on a location distant from the site of the problem. The two residues depicted have nothing in common - they're 80A apart, on different chains, with different identities.

 The phenomenon appears reproducible (for a given starting structure) and is not hardware related - I've tried it on the cluster and on my own local machine (both 64-bit and CUDA); all three are independent NAMD builds. It also seems to have nothing to do with MDFF - the same behaviour occurs whether or not gridforces are on. Each movie shows the initial minimization of the same structure at two different restart points. In one case it's a straightforward restart; in the other, it's after performing an interactive simulation on an unrelated position in the structure and writing a new PDB file. In each case, if allowed to continue the structures settle back to perfectly normal looking configurations by minimisation step 1000-2000, and will go on to uneventful MD. In other cases, similar events have led to trans/cis or chirality flips in local residues.

Everything in the PDB files looks perfectly normal, and there is absolutely nothing visibly untoward in the starting structures. In each case it is clearly the amino acid residue that begins it - subtly for the first 10-30 steps, but quickly mounting into an extreme oscillation along each bond... before settling back completely to normal.

I have no idea how to further work out what is going on here, but I'm at a halt until I know what's causing it and how to avoid it. If it helps, I'm copying below the tcl script used to make the construct using psfgen. The segnames X### and Y### are glycans.

Thanks,

Tristan

# File: dimer_make.tcl
# (2) Embed the psfgen commands in this script
package require psfgen
# Requirements: topology file top_all22_prot.inp in directory toppar
# PDB file complete_xplor.pdb in current directory

# (3) Read topology file
set topdir /home/tristan/ff/CHARMM_c36
topology $topdir/top_all36_prot.rtf
topology $topdir/top_all36_carb.rtf
topology $topdir/stream/toppar_all36_carb_glycopeptide.str
topology $topdir/toppar_water_ions.str
# (4) Build protein segment

foreach i [lsort [glob -directory ./parts -tails *.pdb]] {
  set ii [file rootname $i]
  segment $ii {
    pdb ./parts/$ii.pdb
  }
}

# (5) Apply patches
patch DISU E1:7 E1:26
patch DISU E1:124 E1:152
patch DISU E1:156 E1:179
patch DISU E1:166 E1:185
patch DISU E1:189 E1:198
patch DISU E1:193 E1:204
patch DISU E1:205 E1:213
patch DISU E1:209 E1:222
patch DISU E1:225 E1:234
patch DISU E1:238 E1:250
patch DISU E1:256 E1:277
patch DISU E1:281 E1:295
patch DISU E1:298 E1:302
patch DISU E1:306 E1:327
patch DISU E1:429 E1:462
patch DISU E1:637 E2:853
patch DISU E1:673 E1:676
patch DISU E2:780 E2:789

patch DISU E1:518 F1:518
patch DISU E1:666 F1:666
patch DISU E1:674 F1:674

patch DISU I1:6 I1:48
patch DISU I1:18 I1:61
patch DISU I1:47 I1:52

patch 14bb X25:1 X25:2
patch 16ab X25:1 X25:10
patch 14bb X25:2 X25:3
patch 16ab X25:3 X25:9
patch 13ab X25:3 X25:4
patch 12ba X25:4 X25:5
patch 14bb X25:5 X25:6
patch 14bb X25:4 X25:7
patch 14bb X25:7 X25:8
patch NGLB E1:25 X25:1

patch 14bb X109:1 X109:2
patch 14bb X109:2 X109:3
patch 13ab X109:3 X109:4
patch 12aa X109:4 X109:5
patch 12aa X109:5 X109:6
patch 13ab X109:6 X109:7
patch 16ab X109:3 X109:8
patch 13ab X109:8 X109:9
patch 16ab X109:8 X109:10
patch 12aa X109:10 X109:11
patch NGLB E1:109 X109:1

patch 14bb X218:1 X218:2
patch 14bb X218:2 X218:3
patch 13ab X218:3 X218:4
patch 12aa X218:4 X218:5
patch 16ab X218:3 X218:6
patch 13ab X218:6 X218:7
patch 16ab X218:7 X218:8
patch NGLB E1:218 X218:1

patch 14bb X288:1 X288:2
patch 16ab X288:1 X288:10
patch 14bb X288:2 X288:3
patch 13ab X288:3 X288:4
patch 12ba X288:4 X288:5
patch 14bb X288:5 X288:6
patch 14bb X288:4 X288:7
patch 14bb X288:7 X288:8
patch 16ab X288:3 X288:9
patch NGLB E1:288 X288:1

patch 14bb X391:1 X391:2
patch 14bb X391:2 X391:3
patch 13ab X391:3 X391:4
patch 12aa X391:4 X391:5
patch 16ab X391:3 X391:6
patch 16ab X391:6 X391:7
patch NGLB E1:391 X391:1

patch 14bb X412:1 X412:2
patch 14bb X412:2 X412:3
patch 16ab X412:1 X412:10
patch 13ab X412:3 X412:4
patch 12ba X412:4 X412:5
patch 14bb X412:5 X412:6
patch 14bb X412:4 X412:7
patch 14bb X412:7 X412:8
patch 16ab X412:3 X412:9
patch NGLB E1:412 X412:1

patch 14bb X508:1 X508:2
patch 14bb X508:2 X508:3
patch 13ab X508:3 X508:4
patch 16ab X508:3 X508:5
patch 12aa X508:5 X508:6
patch 12aa X508:6 X508:7
patch NGLB E1:508 X508:1

patch 14bb X581:1 X581:2
patch 14bb X581:2 X581:3
patch 13bb X581:3 X581:4
patch 16ab X581:3 X581:5
patch 12aa X581:5 X581:6
patch 12aa X581:6 X581:7
patch NGLB E1:581 X581:1

patch 14bb X596:1 X596:2
patch 14bb X596:2 X596:3
patch 16ab X596:1 X596:10
patch 13ab X596:3 X596:4
patch 12ba X596:4 X596:5
patch 14bb X596:5 X596:6
patch 14bb X596:4 X596:7
patch 14bb X596:7 X596:8
patch 16ab X596:3 X596:9
patch NGLB E1:596 X596:1

patch 14bb X614:1 X614:2
patch 14bb X614:2 X614:3
patch 16ab X614:1 X614:10
patch 13ab X614:3 X614:4
patch 12ba X614:4 X614:5
patch 14bb X614:5 X614:6
patch 14bb X614:4 X614:7
patch 14bb X614:7 X614:8
patch 16ab X614:3 X614:9
patch NGLB E1:614 X614:1

patch 14bb X874:1 X874:2
patch 14bb X874:2 X874:3
patch 16ab X874:1 X874:10
patch 13ab X874:3 X874:4
patch 12ba X874:4 X874:5
patch 14bb X874:5 X874:6
patch 14bb X874:4 X874:7
patch 14bb X874:7 X874:8
patch 16ab X874:3 X874:9
patch NGLB E2:874 X874:1

patch 14bb X887:1 X887:2
patch 14bb X887:2 X887:3
patch 13ab X887:3 X887:4
patch 12aa X887:4 X887:5
patch 16ab X887:3 X887:6
patch 13ab X887:6 X887:7
patch 16ab X887:7 X887:8
patch NGLB E2:887 X887:1

# (5) Apply patches
patch DISU F1:7 F1:26
patch DISU F1:124 F1:152
patch DISU F1:156 F1:179
patch DISU F1:166 F1:185
patch DISU F1:189 F1:198
patch DISU F1:193 F1:204
patch DISU F1:205 F1:213
patch DISU F1:209 F1:222
patch DISU F1:225 F1:234
patch DISU F1:238 F1:250
patch DISU F1:256 F1:277
patch DISU F1:281 F1:295
patch DISU F1:298 F1:302
patch DISU F1:306 F1:327
patch DISU F1:429 F1:462
patch DISU F1:637 F2:853
patch DISU F1:673 F1:676
patch DISU F2:780 F2:789

patch 14bb Y25:1 Y25:2
patch 16ab Y25:1 Y25:10
patch 14bb Y25:2 Y25:3
patch 16ab Y25:3 Y25:9
patch 13ab Y25:3 Y25:4
patch 12ba Y25:4 Y25:5
patch 14bb Y25:5 Y25:6
patch 14bb Y25:4 Y25:7
patch 14bb Y25:7 Y25:8
patch NGLB F1:25 Y25:1

patch 14bb Y109:1 Y109:2
patch 14bb Y109:2 Y109:3
patch 13ab Y109:3 Y109:4
patch 12aa Y109:4 Y109:5
patch 12aa Y109:5 Y109:6
patch 13ab Y109:6 Y109:7
patch 16ab Y109:3 Y109:8
patch 13ab Y109:8 Y109:9
patch 16ab Y109:8 Y109:10
patch 12aa Y109:10 Y109:11
patch NGLB F1:109 Y109:1

patch 14bb Y218:1 Y218:2
patch 14bb Y218:2 Y218:3
patch 13ab Y218:3 Y218:4
patch 12aa Y218:4 Y218:5
patch 16ab Y218:3 Y218:6
patch 13ab Y218:6 Y218:7
patch 16ab Y218:7 Y218:8
patch NGLB F1:218 Y218:1

patch 14bb Y288:1 Y288:2
patch 16ab Y288:1 Y288:10
patch 14bb Y288:2 Y288:3
patch 13ab Y288:3 Y288:4
patch 12ba Y288:4 Y288:5
patch 14bb Y288:5 Y288:6
patch 14bb Y288:4 Y288:7
patch 14bb Y288:7 Y288:8
patch 16ab Y288:3 Y288:9
patch NGLB F1:288 Y288:1

patch 14bb Y391:1 Y391:2
patch 14bb Y391:2 Y391:3
patch 13ab Y391:3 Y391:4
patch 12aa Y391:4 Y391:5
patch 16ab Y391:3 Y391:6
patch 16ab Y391:6 Y391:7
patch NGLB F1:391 Y391:1

patch 14bb Y412:1 Y412:2
patch 14bb Y412:2 Y412:3
patch 16ab Y412:1 Y412:10
patch 13ab Y412:3 Y412:4
patch 12ba Y412:4 Y412:5
patch 14bb Y412:5 Y412:6
patch 14bb Y412:4 Y412:7
patch 14bb Y412:7 Y412:8
patch 16ab Y412:3 Y412:9
patch NGLB F1:412 Y412:1

patch 14bb Y508:1 Y508:2
patch 14bb Y508:2 Y508:3
patch 13ab Y508:3 Y508:4
patch 16ab Y508:3 Y508:5
patch 12aa Y508:5 Y508:6
patch 12aa Y508:6 Y508:7
patch NGLB F1:508 Y508:1

patch 14bb Y581:1 Y581:2
patch 14bb Y581:2 Y581:3
patch 13bb Y581:3 Y581:4
patch 16ab Y581:3 Y581:5
patch 12aa Y581:5 Y581:6
patch 12aa Y581:6 Y581:7
patch NGLB F1:581 Y581:1

patch 14bb Y596:1 Y596:2
patch 14bb Y596:2 Y596:3
patch 16ab Y596:1 Y596:10
patch 13ab Y596:3 Y596:4
patch 12ba Y596:4 Y596:5
patch 14bb Y596:5 Y596:6
patch 14bb Y596:4 Y596:7
patch 14bb Y596:7 Y596:8
patch 16ab Y596:3 Y596:9
patch NGLB F1:596 Y596:1

patch 14bb Y614:1 Y614:2
patch 14bb Y614:2 Y614:3
patch 16ab Y614:1 Y614:10
patch 13ab Y614:3 Y614:4
patch 12ba Y614:4 Y614:5
patch 14bb Y614:5 Y614:6
patch 14bb Y614:4 Y614:7
patch 14bb Y614:7 Y614:8
patch 16ab Y614:3 Y614:9
patch NGLB F1:614 Y614:1

patch 14bb Y874:1 Y874:2
patch 14bb Y874:2 Y874:3
patch 16ab Y874:1 Y874:10
patch 13ab Y874:3 Y874:4
patch 12ba Y874:4 Y874:5
patch 14bb Y874:5 Y874:6
patch 14bb Y874:4 Y874:7
patch 14bb Y874:7 Y874:8
patch 16ab Y874:3 Y874:9
patch NGLB F2:874 Y874:1

patch 14bb Y887:1 Y887:2
patch 14bb Y887:2 Y887:3
patch 13ab Y887:3 Y887:4
patch 12aa Y887:4 Y887:5
patch 16ab Y887:3 Y887:6
patch 13ab Y887:6 Y887:7
patch 16ab Y887:7 Y887:8
patch NGLB F2:887 Y887:1

foreach i [lsort [glob -directory ./parts -tails *.pdb]] {
  set ii [file rootname $i]
   coordpdb ./parts/$ii.pdb $ii
}
# (9) Guess missing coordinates
guesscoord

# (10) regenerate angles and dihedrals
regenerate angles
regenerate dihedrals

# (10) Write structure and coordinate files
writepsf ionized.psf
writepdb ionized.pdb
# End of psfgen commands

From: Tristan Croll <tristan.croll_at_qut.edu.au<mailto:tristan.croll_at_qut.edu.au>>
Date: Tue, 11 Feb 2014 00:24:37 +0000
To: "namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>" <namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>>
Subject: namd-l: RE: Strange glitches when performing energy minimization in MDFF simulations

I've gone back over the last dozen simulations that involved energy minimization, and I'm seeing similar single-point spikes (ranging from 10^7->10^9 kcal/mol), seemingly randomly distributed in the middle of the minimization process, in about half of them. I can correlate this back to structural deformations in three. In one clear case, I have an arginine residue that looks completely normal at step 1000. Then there's a roughly 10^9 kcal/mol spike in the logfile at step 1062, the arginine is entirely distorted with the water molecules pushed away by about 5A... then by step 5000 the structure is back to normal and MD simulations continue with no sign of trouble. There are no warnings, and no error messages. I've found similar events in two other trajectories, involving entirely unrelated and >50A distant residues. If I have to guess at what's going on, I would say that each of these events correlates with the system losing track of the position of an atom, with the distortions occurring as it snaps back into place.

This is driving me crazy.

From: owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu> [mailto:owner-namd-l_at_ks.uiuc.edu] On Behalf Of Tristan Croll
Sent: Monday, 10 February 2014 4:21 PM
To: namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>
Subject: namd-l: RE: Strange glitches when performing energy minimization in MDFF simulations

Just found this in one of the logfiles. Something's clearly glitching, but where? Is this more likely to be a hardware or a software fault?

ENERGY: 11731207 70944.7907 51205.7902 19764.7817 195.5461 -2487642.8961 216593.6539 0.0000 606.7625 0.0000
  -2128331.5712 0.0000 -2128331.5712 -2128331.5712 0.0000 -7479.5403 -7431.4393 5044132.7704 -7479.5403 -7431.4393
ENERGY: 11731208 72427.9957 51488.6168 19766.0213 214.8399 -2490093.8226 1691328329.2127 0.0000 606.7798 0.0000
  1688982739.6436 0.0000 1688982739.6436 1688982739.6436 0.0000 93196322.9979 93196626.2239 5044132.7704 93196322.9979 93196626.2239
ENERGY: 11731209 70806.8876 51169.4804 19765.2833 192.8327 -2487178.9184 216259.2942 0.0000 606.7590 0.0000
  -2128378.3811 0.0000 -2128378.3811 -2128378.3811 0.0000 -7485.5098 -7444.8419 5044132.7704 -7485.5098 -7444.8419

From: owner-namd-l_at_ks.uiuc.edu<mailto:owner-namd-l_at_ks.uiuc.edu> [mailto:owner-namd-l_at_ks.uiuc.edu] On Behalf Of Tristan Croll
Sent: Monday, 10 February 2014 3:50 PM
To: namd-l_at_ks.uiuc.edu<mailto:namd-l_at_ks.uiuc.edu>
Subject: namd-l: Strange glitches when performing energy minimization in MDFF simulations

Hi,

I'm currently running simulations in which I have a protein complex fitted within a ~10A cryo-EM map, and have been running "pseudo-equilibrium" simulations (with the MDFF constant dialled all the way down to 0.005) to allow it to rattle itself down to an energy minimum with weak guidance from the map. I've found that under these conditions it seems quite safe to remove all secondary structure, cispeptide and chirality extrabonds without anything going awry - except in the specific case where I perform an energy minimization at the start of a continuation run. For some reason, during minimization small, seemingly random portions of the structure (e.g. within a sphere of ~10A diameter) will jump into crazy configurations (bond lengths and angles well away from normal) before settling back to a "standard" arrangement (except sometimes with trans to cis or chirality flips). Apart from the MDFF forces, which shouldn't be nearly strong enough to do this (as I said, everything is very stable in MD without restraints), the rest of the simulation settings are very standard. I'm at a loss to explain what is going on, and this has thrown into question what I thought was >50ns of production run. Can anyone help?

Thanks,

Tristan

NAMD version is 2.9-ibverbs
Using CHARMM-36 parameters
Explicit water, 0.15M sodium chloride
I've checked that the system is neutral

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:22:06 CST