|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
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
|
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
* 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
|
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
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
|
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
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
|
|
|
|
|
Joined: Jun 2004
Posts: 162
Forum Member
|
|
Forum Member
Joined: Jun 2004
Posts: 162 |
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?
|
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
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 ...
|
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
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.
|
|
|
|
|
Joined: Jul 2005
Posts: 18
Forum Member
|
|
Forum Member
Joined: Jul 2005
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.
|
|
|
|
|
Joined: Sep 2003
Posts: 8,660 Likes: 26
Forum Member
|
|
OP
Forum Member
Joined: Sep 2003
Posts: 8,660 Likes: 26 |
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.
|
|
|
|
|
Joined: Nov 2008
Posts: 14
Forum Member
|
|
Forum Member
Joined: Nov 2008
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!
|
|
|
|
|
