CHARMM Development Project Forums CHARMM Community Script Archive bulk dielectric constant from a simulation

Depending on the number of water molecules, 2 ns may not be long enough; a rule of thumb attributed to M. Klein suggests around 500 molecules may be needed for stable MD simulations.

I've recently used the following for 5 ns with a TIP3P water box of 1340 waters and a ca. 34 A edge (set FFX = 32) and get values around 100, which is what others have reported for this water model.

* PME dynamics with 34 A water box
*

stream rtfprm.str  ! RTF, PARAM
stream psfcrd.str  ! PSF, COOR
stream cryst.str   ! CUBIC LATTICE

shake bonh param

open unit 31 write file name dyn.trj
open unit 40 read card name dyn.rea
open unit 41 write card name dyn.res
! NVT ENSEMBLE WITH PRESSURE REPORTING (PMASS 0)
dyna cpt leap restart echeck 250. nstep 500000 nprint 500 -
  iprfrq 5000 atom cdie cutnb 16.0  ctofnb 12.0 -
  wmin 1.5  eps 1.0 inbfrq -1 cutim 16.0 imgfrq -1 -
  ewald pme order 6 spline kappa 0.32 fftx @FFX ffty @FFX fftz @FFX -
  pcons pmass 0.0 pgamma 0.0 pref 1.0  -
  hoover reft 293 first 293 finalt 293 -
  ihtfrq 0 ieqfrq 0 ntrfrq 500 ichecw 0 iasors 1 -
  nsavc 1000 iuncrd 31 iunrea 40 iunwri 41

stop

WARNING: the TCONS keyword was originally included by mistake in the indicated line above; DO NOT use that keyword in conjunction with HOOVER

Thanks! I'll try out the new run parameters and get back to you.

Have you tried the new velocity Verlet integrator (VV2)? The last time I simulated a water box with it, I got a large persistent net dipole that wouldn't go away, and very small fluctuations.

I have not used VV2

Hey there,

I just wanted to say thanks for your tips on calculating the dielectric constant. The Ewald summation was what did the trick; without it, the dielectric constant was too low (and rightfully so - the Coulomb interaction between molecules had been weakened). I got a dielectric constant of 71.8 after a 10ns simulation of 216 SPC/E water molecules. I also got a dielectric constant 54 of TIP4P. Now I can move onto developing polarizable models.

Thanks again!

Hi Rick,

Could you please tell me the expression from the dielectric constant to the dipole values of each frame in a trajectory?

Thanks.

---
Jeffrey

I'm away from the office, but I believe this was a key reference--

R. D. Mountain, D. Thirumalai, "Ergodic Measures for the Simulation of Dielectric Properties of Water", Comp. Phys. Comm. 62, 352-359 (1991).

Hi all,

I am new to CHARMM. I am trying to get a vector time series of the system dipole.However, no useful info is written in the boxdipol.dat

Here is the input:

READ .RTF AND .PRM
READ SEQUENCE AND COORDINATES
LONE PAIR FACILITY
SET PBC USING CRYSTAL

! box dipole
open read unit 32 unformatted name ../tip4p-nvt/output/dyna.dcd

traj quer unit 32
traj first 32 nunit 1 skip ?SKIP
open write unit 33 card name boxdipol.dat
echu 33

set i 1
label loop
traj read

coor dipole sele all end

write title unit 33
*

incr i by 1
if i le ?NFILE goto loop
!End the process
close unit 33

stop

It is my understanding from the output that the trajectory was read correctly and dipole moments are calculated but the information is not being written in boxdipol.dat

Here are some parts from the output:

CHARMM> traj quer unit 32

READING TRAJECTORY FROM UNIT 32
NUMBER OF COORDINATE SETS IN FILE: 500
NUMBER OF PREVIOUS DYNAMICS STEPS: 301000
FREQUENCY FOR SAVING COORDINATES: 1000
NUMBER OF STEPS FOR CREATION RUN: 500000

NUMBER OF DEGREES OF FREEDOM: 1533
NUMBER OF FIXED ATOMS: 0
THE INTEGRATION TIME STEP (PS): 0.0010
THE FILE CONTAINS CRYSTAL UNIT CELL DATA
THE FILE DOES NOT CONTAIN 4-D DATA

CHARMM> traj first 32 nunit 1 skip ?SKIP
RDCMND substituted energy or value "?SKIP" to "1000"
TRAJ: INITIATING READ OF A TRAJECTORY, OPTIONS;
FIRSTU = 32 NUNIT = 1 SKIP = 1000

CHARMM> open write unit 33 card name boxdipol.dat
VOPEN> Attempting to open::boxdipol.dat::
OPNLGU> Unit 33 opened for WRITE access to /waterbox/analysis/boxdipol.dat

CHARMM> echu 33

CHARMM> set i 1
Parameter: I <- "1"

CHARMM> label loop

CHARMM> traj read

READING TRAJECTORY FROM UNIT 32
NUMBER OF COORDINATE SETS IN FILE: 500
NUMBER OF PREVIOUS DYNAMICS STEPS: 301000
FREQUENCY FOR SAVING COORDINATES: 1000
NUMBER OF STEPS FOR CREATION RUN: 500000

TITLE> * A BOX OF 256 TIP4P WATER MOLECULES
TITLE> * 500 PS NVT DYNAMICS - PRODUCTION
TITLE> *
***** WARNING ***** BEGIN= 0 Was not specified. It has been set to: 301000

CHARMM> coor dipole sele all end
SELRPN> 1024 atoms have been selected out of 1024
THE TOTAL CHARGE OF SELECTED ATOMS IS: 0.000000
DIPOLE MOMENT ABOUT CENTER OF GEOMETRY (DEBYES) : -95.270272 13.593494 -36.363997 102.876373

CHARMM> write title unit 33
RDTITL> *
RDTITL> No title read.

CHARMM> incr i by 1
Parameter: I <- "2"

CHARMM> if i le ?NFILE goto loop
RDCMND substituted energy or value "?NFILE" to "500"
Comparing "2" and "500".
IF test evaluated as true. Performing command
CHARMM> traj read

CHARMM> coor dipole sele all end
SELRPN> 1024 atoms have been selected out of 1024
THE TOTAL CHARGE OF SELECTED ATOMS IS: 0.000000
DIPOLE MOMENT ABOUT CENTER OF GEOMETRY (DEBYES) : -111.358020 22.111115 -51.917974 124.839842

Thanks in advance for your help

You have connected ("opened") unit/handle 33 to the file boxdipol.dat, but there are no commnands in your input file that actually write anything to this unit.

How about using the ECHO command (for which the outputfile was defined with the ECHU command) instead of the WRITE? It would also help if you actually write something to the file ... eg:
ECHO hi this is a dipole moment: ?XDIP ?YDIP ?ZDIP

(see subst.doc for CHARMM internal variables)

or use the CORREL module

Great, thanks!

How can I get a histogram of dipole moments for example on x dimension (position)? In other words, how can get the corresponding positions of molecules after calculating their dipoles?