You are here

temperature rises during MD run

Forums: 

Hello, I am a new Unres user, and I try to simulate the cohesin protein (pdb code 1AOH) in UNRES (e0ll2y force field was used). I minimized the protein using the following input file:

1AOH minimalizacja energii start ze struktury udokladnionej :
SEED=-1248960 minimize cart pdbref refstr rescale_mode=2 out1file pdbstart &
pdbout
WSC=1.00000 WSCP=1.23315 WELEC=0.84476 WBOND=1.00000 WANG=0.62954 &
WSCLOC=0.10554 WTOR=1.84316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000 &
WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000 &
WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.25 &
CUTOFF=7.00000 WCORR4=0.00000
1AOH.pdb
0
0

The minimisation was successful (RMSD about 2 A from PDB structure), so I used the minimized PDB to run MD simulations with the following input file:

1AOH dynamika Langevina
SEED=-3059743 PDBREF refstr t_bath=300 tbf unres_pdb dt=0.001 &
nstep=1000000 ntwe=100 ntwx=1000 lang=1 pdbout mdpdb pdbstart

WSC=1.00000 WSCP=1.23315 WELEC=0.84476 WBOND=1.00000 WANG=0.62954 &
WSCLOC=0.10554 WTOR=1.84316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000 &
WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000 &
WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.25 &
CUTOFF=7.00000 WCORR4=0.00000
1AOH_min_GB000.pdb
0
0

During the simulation temperature and radius of gyration rise to very high values (like few thousand degrees). What is the possible cause of that? How can I prevent that behaviour?

Start from UNRES_PDB can give weird results. In the new version of UNRES, we introduced the PREMINIM option which minimizes the energy before MD. However, this feature is still absent in the canonical PLGRID version and I'd recommend the following procedure:

1. Do a short (10,000 step) run with very small time step (e.g., dt=0.01 or even dt=0.001). Specify the RESTART keyword on the MD card.

2. Watch the potential energy, if it drops below 0 or at least to the order of 10**2 - 10**3 you can probably restart the calculaltion with regular timestep (dt=0.1).

3. If the potential energy is still high, you need more than one such preparatory run. Don't be confused that the kinetic temperature gets high; it will be brought to normal after you restart the run (also, you can specify reset_vel=1000 or so to have the velocities reset to random numbers corresponding to the thermostat temperature every 1000 timesteps or so).

The above procedure works like minimization or even better.

If your molecule still blows up, you may need to impose distance restraints on part of the structure; the restraints are calculated from the experimental PDB structure provided (section 8.1.13 of the manual provided at www.unres.pl/unres)

I noticed that there was a bug in the energy term weights in the examples on PLGRID; if you use the enhanced 1le11l2y force field with wsccorr=0.25, then you need to use WTOR=1.34316 and not WTOR=1.84316. Please correct.

Adam

Thank you for the quick answer, the problem was actually pretty simple: In the input file I did not make a separate line for control data and MD data, but used the same line for both, which produced strange results. In the correct simulation temperature is also high at the very beginning (because, as suggested, I did not use minimized structure), but quickly goes down to 300 K. The correct input file (for the classic, not enhanced e0ll2y force field) is:

1AOH dynamika Langevina (klasyczne pole e0ll2y)
SEED=-3059743 PDBREF MD pdbstart pdbout
nstep=1000000 ntwe=100 ntwx=10000 dt=0.01 damax=10000.0 lang=1 scal_fric=0.01 &
t_bath=300 mdpdb
WSC=1.00000 WSCP=1.23315 WELEC=0.84476 WBOND=1.00000 WANG=0.62954 &
WSCLOC=0.10554 WTOR=1.84316 WTORD=1.26571 WCORRH=0.19212 WCORR5=0.00000 &
WCORR6=0.00000 WEL_LOC=0.37357 WTURN3=1.40323 WTURN4=0.64673 WTURN6=0.00000 &
WVDWPP=0.23173 WHPB=1.00000 WSCCOR=0.0 CUTOFF=7.00000 WCORR4=0.00000
1AOH.pdb
0
0

Best regards,
Lukasz