There was a problem in my input file that led to spurious results. It wasn't any problem with the CHARMM program.
But I wanted to point out the following.
This is what the CHARMM tutorial says about NVT dynamics:
"PCONS: Turns on pressure reporting. The pressure of the system will be held constant unless PMASS is set to 0
(in this case, pressure statistics are still displayed, but constant pressure is not enforced and an NVT ensemble is obtained)."
This is a correct, but albeit somewhat facile statement.
Below is a snippet of a CHARMM script for a solvated decaalanine solvated in a rhombic dodecahedron simulation cell. This was part of a tutorial for a QM and MD workshop presented in May. The decaalanine was constructed as an idealized alpha-helix and obtained from a CHARMM script in one of the earlier versions of CHARMM. The solvation was performed as given in the Full Example.
! NVT Dynamics
dyna leap cpt restart - ! restart the run from our heating restart, use CPT dynamics
timestep 0.001 nstep 1000000 nprint 1000 - ! do 1.0 ns in 1 fs time step
iunrea 31 iunwri 41 iuncrd 51 isvfrq 1000 nsavc 1000 - ! units for reading/writing restarts and coordinate trajectory
hoover reft 300.0 tmass @tmass - ! Hoover thermostat, constant temp - 300 K, tmass as above
pcons pref 1.0 pmass 0.5 pgamma 0 - ! NVT simulation, pmass 0.5 (need pmass~0), collision frequency 0
ihtfrq 0 ieqfrq 0 ntrfrq 1000 - ! temperature under Hoover's control, cancel rotation/translation every 1000 steps
iasors 1 iasvel 1 iscvel 0 ichecw 0 - ! since we're restarting and ihtfrq = ieqfrq = 0, these don't really matter
echeck 100. ! check energy diff as before
The above partial script is similar to that provided in the tutorial, but modified for use in this simulation. Notice that pmass is NOT equal to zero for a perfect NVT simulation. When pmass=0, the system becomes extremely unstable. To approach the value of pmass=0.50 required four separate 1 ns simulations:
1. Step 1
! calculate masses for the pressure and temperature pistons (pmass & tmass)
scalar wmain = mass
scalar mass stat select segid HLX end
calc pmass = ?stot * 0.02
calc tmass = ?stot * 0.20
calc pmass = int ( ?stot / 50.0 )
calc tmass = @pmass * 10
2. Step 2, Restart from Step 1
set pmass = 10.0
3. Step 3, Restart from Step 2
set pmass = 1.0
4. Step 4, Restart from Step 3
set pmass = 0.5
With pmass = 0.5, the cell parameters/volume are quite stable for 1 ns, but for a larger system and longer simulation time (~ 50 ns) instabilities may arise.
I apologize for initially thinking that CHARMM has a problem, when it was my input error in the image command.