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
Page 1 of 3 1 2 3 >

Moderator:  chmgr, John Legato, petrella