Page 1 of 3 1 2 3 >
Topic Options
#4290 - 11/17/04 09:16 PM bulk dielectric constant from a simulation
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
The following CHARMM script and Fortran program are used together to compute the dielectric constant of cubic systems; I've been using them to evaluate developmental parameters on simulations of neat fluids. The script is used to create a vector time series of the system dipole using COOR DIPOLE (corman.doc) in a loop over frames. The Fortran program computes eps from that time series, and outputs the cumulative value of eps vs. time as a text data file suitable for plotting. There's some unit conversion going on the program, from AKMA into SI units.
_________________________
Rick Venable
computational chemist


Top
#4291 - 11/17/04 09:29 PM Re: bulk dielectric constant from a simulation [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
* box dipole time series
*

! READ RTF, PARAM FILE
stream rtfprm.str
! READ PSF AND INITIAL (T=0) COORDS
stream psfcrd.str
! SETUP CUBIC WITH CRYST, TO GET XTLA FROM TRJ FILES
stream initcrys.str

! NO. OF FILES
set nf 20
! NO. OF TOTAL FRAMES
calc nt = @NF * 100

! SETUP TRJ FILES FOR READING IN A LOOP
set k 1
label oplp
calc u 20 + @K
open unit @U write file name dyn@K.trj
incr k by 1
if k .le. @NF goto oplp
traj first 21 nunit @NF

! OPEN DATA FILE, INIT LOOP COUNTER
open unit 12 write card name bxdipol.dat
set j 1
label jlp
! READ A FRAME, COMPUTE DIPOLE OF ALL MOLECULES
traj read
coor dipole sele all end
! WRITE TIME, DIPOLE VECTOR, CELL EDGE, DIPOLE MAGNITUDE
write title unit 12
* @J. ?XDIP ?YDIP ?ZDIP ?XTLA ?RDIP
*
! INCREMENT COUNTER, LOOP EXIT TEST
incr j by 1
if j le @NT goto jlp

close unit 12
stop

Top
#4292 - 11/17/04 09:41 PM Re: bulk dielectric constant from a simulation [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
      program dielec
C COMPUTE DIELECTRIC FROM TIME SERIES OF BOX DIPOLE MAGNITUDE
C use Coulombs, meters, Joules
integer i,j,k, nt,mxt
character*80 fn
parameter (mxt=10000)
real*8 dp(mxt), d2(mxt), ti(mxt),ed(mxt), bx(mxt),by(mxt),bz(mxt)
real*8 adp,ad2, pi,tem, tv,av, mx,my,mz, rdi,rd2, sd2

C KB == Boltzmann, DC == Debye to C.m, EPS0 == vac permittivity
REAL*8 KB, DC, EPS0
PARAMETER (KB=1.3805d-23, DC=3.336d-30, EPS0=8.85419d-12)

pi = acos(-1.d0)
tem = 293.15d0

k=iargc()
if (k.ne.2) then
write(6,*) 'dielec.ix DIPOLE.DAT CUMDIEL.DAT'
write(6,*) ' DIPOLE.DAT input; time x y z edge dipole_mag'
write(6,*) ' CUMDIEL.DAT output; time cum_dielec'
stop 'and try again'
else
call getarg(1,fn)
open(1,file=fn,form='formatted')
call getarg(2,fn)
open(2,file=fn,form='formatted')
endif

nt = 1
10 read(1,*,end=20,err=20) ti(nt),bx(nt),by(nt),bz(nt),ed(nt),dp(nt)
nt = nt + 1
if (nt.le.mxt) goto 10
20 continue
if (nt.lt.mxt) nt = nt-1
fac = (dc*dc) / (3.d0*kb*tem*eps0)
write (6,*) 'Scale factor = ', fac

tv = 0.d0
mx = 0.d0
my = 0.d0
mz = 0.d0
sd2 = 0.d0
do k=1,nt

mx = mx + bx(k)
my = my + by(k)
mz = mz + bz(k)
rdi = sqrt(mx**2 + my**2 + mz**2)
sd2 = sd2 + bx(k)**2 + by(k)**2 + bz(k)**2
tv = tv + (1.d-10*ed(k))**3

av = tv / dble(k)
adp = rdi / dble(k)
ad2 = sd2 / dble(k)
write (2,'(2F11.1)') ti(k),1.d0 + fac*(ad2-adp**2)/av

enddo

write (6,*) 'Final eps=', 1.d0 + fac*(ad2-adp**2)/av,
1 ' < M**2 > =',ad2, ' < M > =',adp

stop
end

Top
#4293 - 04/23/07 11:08 AM Re: bulk dielectric constant from a simulation [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
There's a bug in the above CHARMM script; the trajectory files are incorrectly opened for write access. Change

open unit @U write file name dyn@K.trj

to

open unit @U read file name dyn@K.trj
_________________________
Rick Venable
computational chemist


Top
#4294 - 08/17/07 03:18 PM Re: bulk dielectric constant from a simulation [Re: rmv]
noskov Offline
Forum Member

Registered: 06/18/04
Posts: 162
Loc: Calgary, AB
Rick,

First, let me thank you for a scripts. Very useful, indeed. However, it seems that in versions of CHARMM after c32 the output of correl series and dipole command for X,Y and Z vectors was somehow adjusted to give time series in Debye and not eA as it used to be. so conversion factor for more recent charmm versions would be bit different and factor 4.8 needs to be presented.
Am I right?
_________________________
Sergei Noskov snoskov@ucalgary.ca Dept.BioSciences University of Calgary

Top
#4295 - 08/17/07 06:27 PM Re: bulk dielectric constant from a simulation [Re: noskov]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
Ugh. I'll have to look at the code, I guess. Dealing with electrostatic constants is enough of a pain, without arbitrary changes in the units ...

Top
#4296 - 10/21/08 06:26 PM Re: bulk dielectric constant from a simulation [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
I've used this code more recently and gotten reasonable answers; my comments in the above Fortran code seem to indicate that the input data is expected to be in Debye units.

Top
#4297 - 11/12/08 05:40 PM Re: bulk dielectric constant from a simulation [Re: rmv]
Mandy Offline
Forum Member

Registered: 07/15/05
Posts: 18
can I use the Fortran program to calculate dielectric constant for protein-water system? where the time series of the system dipole still be generated using COOR DIPOLE, from a trajectory file.
thanks.

Top
#4298 - 11/12/08 06:23 PM Re: bulk dielectric constant from a simulation [Re: Mandy]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
I believe it should work, as long as the system does not have a net charge. In theory, the DIPOLE time series could also be used, but that feature was added after I write the above script and program.

Top
#4299 - 11/13/08 04:17 PM Re: bulk dielectric constant from a simulation [Re: rmv]
leeping Offline
Forum Member

Registered: 11/13/08
Posts: 14
I've been trying for the past week or so to reproduce the dielectric constant for water. I'm using the SPC/E model as implemented in CHARMM. However, I am consistently getting values for the dielectric constant that are too low (~55 instead of the reported 70). My dielectric constants are calculated from rather long trajectories (~2ns) and the values appear to be converged.

In your experience, which run parameters have the greatest influence on calculating the dielectric constant? Which run parameters are insignificant?

If you happen to have a set of input files that correctly reproduces the dielectric constant for water (any model will do), may I get it from you?

Thanks a lot!

Top
#4300 - 11/13/08 05:50 PM Re: bulk dielectric constant from a simulation [Re: leeping]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
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


Edited by rmv (02/28/09 03:07 PM)

Top
#4301 - 11/13/08 10:09 PM Re: bulk dielectric constant from a simulation [Re: rmv]
leeping Offline
Forum Member

Registered: 11/13/08
Posts: 14
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.

Top
#4302 - 11/13/08 11:20 PM Re: bulk dielectric constant from a simulation [Re: leeping]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
I have not used VV2

Top
#4303 - 11/19/08 11:25 PM Re: bulk dielectric constant from a simulation [Re: rmv]
leeping Offline
Forum Member

Registered: 11/13/08
Posts: 14
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!

Top
#4304 - 02/28/09 07:21 AM Re: bulk dielectric constant from a simulation [Re: rmv]
jeffrey Offline
Forum Member

Registered: 03/18/07
Posts: 148
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

Top
#4305 - 02/28/09 01:41 PM Re: bulk dielectric constant from a simulation [Re: jeffrey]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
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).

Top
#21244 - 06/13/09 02:31 PM Re: bulk dielectric constant from a simulation [Re: rmv]
TGK Offline
Forum Member

Registered: 05/22/09
Posts: 17
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

Top
#21245 - 06/13/09 02:52 PM Re: bulk dielectric constant from a simulation [Re: TGK]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4741
Loc: ~ 59N, 15E
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
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#21246 - 06/13/09 04:11 PM Re: bulk dielectric constant from a simulation [Re: lennart]
TGK Offline
Forum Member

Registered: 05/22/09
Posts: 17
Great, thanks!

Top
#22444 - 10/09/09 02:18 AM Re: bulk dielectric constant from a simulation [Re: TGK]
TGK Offline
Forum Member

Registered: 05/22/09
Posts: 17
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?


Edited by TGK (10/09/09 02:42 AM)

Top
#22457 - 10/09/09 01:47 PM Re: bulk dielectric constant from a simulation [Re: TGK]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
Within CORREL (correl.doc), you can create matching time series of a molecule's dipole moment vector or x-component (DIPOLE time series) and it's center of mass vector or x-component (ATOM time series with MASS keyword and an atom selection).

The rest would require a user-written program or procedure, but it would be a fairly simple one.
_________________________
Rick Venable
computational chemist


Top
Page 1 of 3 1 2 3 >

Moderator:  chmgr, John Legato, petrella