Previous Thread
Next Thread
Print Thread
Page 1 of 3 1 2 3
Joined: Sep 2003
Posts: 8,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
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,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
* 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,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
      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,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
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
Offline
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?


Sergei Noskov snoskov@ucalgary.ca Dept.BioSciences University of Calgary
Joined: Sep 2003
Posts: 8,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
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,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
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
M
Forum Member
Offline
Forum Member
M
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,605
Likes: 24
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
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
L
Forum Member
Offline
Forum Member
L
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!

Page 1 of 3 1 2 3

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.010s Queries: 35 (0.006s) Memory: 0.7841 MB (Peak: 0.8790 MB) Data Comp: Off Server Time: 2022-09-28 00:20:44 UTC
Valid HTML 5 and Valid CSS