<html>

The problem with Ewald sums is that, although it has other great virtues,

it is not strictly conservative, and with the approximations of PME that

becomes a serious problem and pretty well precludes the use of NVE

boundary conditions (at least that is our experience).&nbsp; We

considered some of this using just pure water as a test system

in:<br><br>

<font face="MS Serif, Geneva">Chiu, S.-W., Clark, M., Subramaniam, S.,

and E. Jakobsson. 2000. Collective motion artefacts arising in

long-duration molecular dynamics simulations.&nbsp;&nbsp; <i>J. Comp.

Chem. </i>21:121-131<br><br>

</font>Eric <br><br>

At 01:49 PM 3/4/2003 -0500, you wrote:<br>

<blockquote type=cite class=cite cite>Hi, I am trying to simulate a

protein in flexible spc water and I believe I have very conservative

parameters. I however see a drift in the total energy, and although it's

not too large it will become large as I need the simulation to be several

ns long. If you plot the total energy you see that the system is

consistently heating up. I have equilibrated the system for more than 100

ps doing npt and before that using position restraints,<br>

so everything should be fine ...<br><br>

time step is 1 fs ,&nbsp; 1nm for rvdw, fourierspacing is 0.12 for pme ,

i update the neighbor list every step nstlist =1, etc ...<br>

the mdp file is attached.<br>

Do you think it's the PME? if so,&nbsp; what should i use for fourier

spacing and tolerance to make it not drift? I am running single precision

though.<br>

Thanks!!!<br>

Claudio.<br><br>

here is the output from g_energy statistics over 380 ps<br><br>

Energy&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

Average&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; RMSD&nbsp;&nbsp;&nbsp;&nbsp;

Fluct.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Drift Tot-Drift<br>

-------------------------------------------------------------------------------<br>

Bond&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

15717.5&nbsp;&nbsp;&nbsp; 254.829&nbsp;&nbsp;&nbsp; 252.943&nbsp;

-0.282489 -107.205<br>

Angle&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

10366&nbsp;&nbsp;&nbsp; 194.383&nbsp;&nbsp;&nbsp; 193.899&nbsp;&nbsp;

0.125137&nbsp; 47.4897<br>

Proper

Dih.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

122.033&nbsp;&nbsp;&nbsp; 12.9951&nbsp;&nbsp;&nbsp; 12.9539

0.00944353&nbsp;&nbsp; 3.58383<br>

Ryckaert-Bell.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

1424.82&nbsp;&nbsp;&nbsp; 52.2058&nbsp;&nbsp;&nbsp; 52.0629 -0.0352289

-13.3694<br>

LJ-14&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

1983.58&nbsp;&nbsp;&nbsp; 41.7575&nbsp;&nbsp;&nbsp; 41.6305

-0.0297046&nbsp; -11.2729<br>

Coulomb-14&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

11277.9&nbsp;&nbsp;&nbsp; 55.3076&nbsp;&nbsp;&nbsp; 53.6154&nbsp;&nbsp;

0.123927&nbsp; 47.0305<br>

LJ

(SR)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

23216.1&nbsp;&nbsp;&nbsp; 415.996&nbsp;&nbsp;&nbsp; 415.932

-0.0663132&nbsp; -25.1659<br>

Disper.

corr.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

-776.671&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

0&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

0&nbsp;&nbsp;&nbsp; 0<br>

Coulomb

(SR)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

-161371&nbsp;&nbsp;&nbsp; 627.813&nbsp;&nbsp;&nbsp; 626.229&nbsp;&nbsp;

0.406744&nbsp;&nbsp; 154.36<br>

Coulomb

(LR)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

-36995&nbsp;&nbsp;&nbsp; 51.9416&nbsp;&nbsp;&nbsp; 50.1115&nbsp;

-0.124746&nbsp;&nbsp; -47.3413<br>

Potential&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

-135035&nbsp;&nbsp;&nbsp; 251.587&nbsp;&nbsp;&nbsp; 251.203&nbsp;&nbsp;

0.126766&nbsp;&nbsp;&nbsp; 48.1079<br>

Kinetic

En.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

37008.7&nbsp;&nbsp;&nbsp; 245.956&nbsp;&nbsp;&nbsp; 245.244&nbsp;&nbsp;

0.170733&nbsp;&nbsp;&nbsp; 64.7933<br>

Total

Energy&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

-98026.1&nbsp;&nbsp;&nbsp; 35.0014&nbsp;&nbsp;&nbsp; 12.7625&nbsp;&nbsp;

0.297498&nbsp;&nbsp;&nbsp; 112.901<br><br>

<br>

-- <br>

-----------------------------------------------<br>

Claudio Javier Margulis,

Ph.D.&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

<br>

Chemistry

Department&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

<br>

Columbia

University&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

<br>

Web

<a href="http://www.chem.columbia.edu/~claudiom/%A0%A0%A0%A0%A0%A0%A0" eudora="autourl">http://www.chem.columbia.edu/~claudiom/&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

</a> <br>

Work phone (212)

854-8470&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

<br>

-----------------------------------------------<br><br>

<br><br>

;;here no constraints and<br>

; nslist=1 nstcomm removed as parameter<br>

;notice that by default gromacs makes center of mass velocity

removal<br>

; here we are trying flexible spc to see if we get better energy

consevation<br>

;<x-tab>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</x-tab>User spoel

(236)<br>

;<x-tab>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</x-tab>Wed Nov&nbsp; 3

17:12:44 1993<br>

;<x-tab>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</x-tab>Input 

file<br>

;<br>

title&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; Yo<br>

cpp&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; /lib/cpp<br>

define&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; -DFLEX_SPC<br>

constraints&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =&nbsp;

none<br>

integrator&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =&nbsp;

md<br>

dt&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 0.001<x-tab>&nbsp;&nbsp;&nbsp;&nbsp;</x-tab>; ps !<br>

nsteps&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 50000000<x-tab>&nbsp;</x-tab>; 50 ns<br>

;nstcomm&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 1 should not be there?<br>

nstxout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 50000<br>

nstxtcout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 1500<br>

nstvout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 50000<br>

nstfout&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 0<br>

nstlog&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 5000<br>

nstenergy&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 1500<br>

nstlist&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 1<br>

ns_type&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; grid ;n-list parameter<br>

rlist&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 1.0 ;cutoff for the short range neighbor list<br><br>

<br>

; OPTIONS FOR ELECTROSTATICS AND VDW = <br>

; Method for doing electrostatics = <br>

coulombtype&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= PME<br>

; Method for doing Van der Waals = <br>

vdw-type&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= Cut-off<br>

; cut-off lengths&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = <br>

rvdw-switch&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 0<br>

rvdw&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 1.0<br>

; Apply long range dispersion corrections for Energy and Pressure = 

<br>

DispCorr&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= EnerPres<br>

; Spacing for the PME/PPPM FFT grid = <br>

fourierspacing&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 0.12<br>

; FFT grid size, when a value is 0 fourierspacing will be used = <br>

fourier_nx&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 0<br>

fourier_ny&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 0<br>

fourier_nz&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 0<br>

; EWALD/PME/PPPM parameters = <br>

pme_order&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 4<br>

ewald_rtol&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 1e-05<br>

ewald_geometry&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= 3d<br>

epsilon_surface&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; =

0<br>

optimize_fft&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= yes<br>

; Nose-Hoover temperature coupling is off<br>

Tcoupl&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= no<br>

tc-grps<x-tab>&nbsp;</x-tab>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=<br>

tau_t&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=<br>

ref_t&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=<br>

&nbsp;Energy monitoring<br>

energygrps<x-tab>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</x-tab>&nbsp;&nbsp;&nbsp;

=&nbsp; Protein SOL Na<br>

; Pressure coupling is off<br>

Pcoupl&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= no<br>

Pcoupltype&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = <br>

tau_p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= <br>

compressibility&nbsp;&nbsp;&nbsp;&nbsp; = <br>

ref_p&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

= <br>

; Generate velocites is on at 300 K.<br>

gen_vel&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; yes<br>

gen_temp&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 300.0<br>

gen_seed&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;

=&nbsp; 173529</blockquote>

<x-sigsep><p></x-sigsep>

---------------------------------<br>

Eric Jakobsson, Ph.D.<br>

Professor, Department of Molecular and Integrative Physiology, and of

Biochemistry<br>

Senior Research Scientist, National Center for Supercomputing

Applications<br>

Professor, Beckman Institute for Advanced Science and Technology<br>

4021 Beckman Institute, mc251<br>

405 N. Mathews Avenue<br>

University of Illinois, Urbana, IL 61801<br>

ph. 217-244-2896&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; fax

217-244-2909<br><br>

<br>

</html>