PROGRAM B2
IMPLICIT NONE
! This program calculates the second virial coefficient of two
! groups of molecules.
! Define inputs: Nc,Np = number of beads in group 1 and 2
! Xc,Yc,Zc = coordinates of beads in group1.
! Xp,Yp,Zp = coordinates of beads in group2.
! XcRot,YcRot,ZcRot = Rotated coordinates of
! beads in group1.
! XpRot,YpRot,ZpRot = Rotated coordinates of
! beads in group2.
! Mc, Mp = Mass of beads in group 1 and 2
! Energy = energy between groups.
! FFactor = exp(U/BT)-1
! TotalF = sum of Ffactors.
! AveFfactor = average of Ffactors
! EPS, SIG = epsilon and sigma matrices.
! XCcm,YCcm,ZCcm,XPcm,YPcm,ZPcm = center of
! mass coordinates for group c and p.
! Kappa, Lambda = correction to bertelot and
! lorentze rule for epsilon and sigma.
! ThetaC,ThetaP = angle of rotation for group c,p
! r, maxR = distance,max.distance between groups.
Logical first_inc ! records the first increase
! in f factor.
integer i,j, h, k, Nc, Np, seed1, seed2, seed3, seed4
integer axisC, axisP
integer numtrials ! number of rotation trials.
integer numrot1,numrot2,numrot3,numrot4 ! Number of rotations for each segment.
real disrot1,disrot2,disrot3,disrot4 ! distance to start new number of rotation.
character (Len = 20) fileinput,fileresult
character (Len = 20) Name ! name of type of molecules.
real, dimension(:), allocatable :: Xc(:), Yc(:), Zc(:)
real, dimension(:), allocatable :: XcRot(:), YcRot(:), ZcRot(:)
real, dimension(:), allocatable :: XpRot(:), YpRot(:), ZpRot(:)
real, dimension(:), allocatable :: Xcnew(:), Ycnew(:), Zcnew(:)
real, dimension(:), allocatable :: Xp(:), Yp(:), Zp(:)
real, dimension(:), allocatable :: Xpnew(:), Ypnew(:), Zpnew(:)
real, dimension(:), allocatable :: Mp(:), Mc(:)
real, dimension(:), allocatable :: Energy(:)
real, dimension(:), allocatable :: Ffactor(:)
real, dimension(:), allocatable :: TotalF(:)
real, dimension(:), allocatable :: AveFfactor(:)
real, dimension(:,:,:), allocatable :: EPS(:,:,:), SIG(:,:,:)
real XPcm, YPcm, ZPcm, XCcm, YCcm, ZCcm ! cm coordinates.
real OldXPcm, OldYPcm, OldZPcm
real, dimension(:,:), allocatable :: Kappa(:,:), lambda(:,:)
double precision :: thetaC, thetaP
real, parameter :: pi=3.141592654
real maxR ! maximum distance between groups.
real r,r2 ! r and r square.
real Norm_inc ! Length of normal increments.
real increment, dt ! distance increment.
real small_inc
real change_R ! changes from small to large increments.
real change_rot ! changes number of rotation trials.
real ran2
real T ! temperature (K)
real MidR
real LastR ! Longest distance actually used in calculation.
real diff ! difference between Fold and Fnew.
real IntVal ! value of integral.
real b2val
real b2final
integer kappcount ! counter for kappa
integer, dimension(:), allocatable :: TYPEp(:), TYPEc(:)
integer Nham ! number of hailtonians
integer LJgroups ! number of LJ groups.
integer num_ij ! number of different kappa's,lambda's.
integer num_points ! number of points in output file.
seed1= -123
seed2= -345
seed3= -124
seed4= -653
num_points = 0
Ffactor = 0.0
TotalF = 0.0
kappcount = 1
first_inc =.true.
fileinput = 'rvsf.dat'
fileresult = 'b2result.dat'
!******************* Read Data ***********************
open(10, file = 'b2data.dat')
read(10,*)
read(10,*)
read(10,*) Name
read (10,*)
read (10,*) T
read(10,*)
read(10,*) Nc
Allocate(Xc(Nc));Allocate(Yc(Nc));Allocate(Zc(Nc));Allocate(Mc(Nc))
Allocate(Xcnew(Nc));Allocate(Ycnew(Nc));Allocate(Zcnew(Nc))
Allocate(XcRot(Nc));Allocate(YcRot(Nc));Allocate(ZcRot(Nc))
read (10,*)
read (10,*) Np
Allocate(Xp(Np));Allocate(Yp(Np));Allocate(Zp(Np));Allocate(Mp(Np))
Allocate(Xpnew(Np));Allocate(Ypnew(Np));Allocate(Zpnew(Np))
Allocate(TYPEc(Nc));Allocate(TYPEp(Np))
Allocate(XpRot(Np));Allocate(YpRot(Np));Allocate(ZpRot(Np))
!Allocate(XpUnRot(Np));Allocate(YpUnRot(Np));Allocate(ZpUnRot(Np))
read (10,*)
read (10,*)
do i=1,Nc
read(10,*) typec(i), Mc(i), xc(i), yc(i), zc(i)
end do
read (10,*)
read (10,*)
do i=1,Np
read(10,*) typep(i), Mp(i), xp(i), yp(i), zp(i)
end do
read (10,*)
read (10,*)
read (10,*) disrot1, numrot1
read (10,*) disrot2, numrot2
read (10,*) disrot3, numrot3
read (10,*) disrot4, numrot4
read (10,*)
read (10,*) Nham
Allocate(Energy(Nham)) ! allocate energy for each hamiltonian.
Allocate(TotalF(Nham))
Allocate(Ffactor(Nham))
!Allocate(Fnew(Nham))
!Allocate(Fold(Nham))
Allocate(AveFfactor(Nham))
read (10,*)
read (10,*) LJgroups
Allocate(EPS(LJgroups,LJgroups,Nham)) !allocate epsilon.
Allocate(SIG(LJgroups,LJgroups,Nham)) !allocate sigma.
read (10,*)
read (10,*)
read (10,*)
do k=1,Nham
read(10,*) (sig(i,i,k), i=1,LJgroups)
end do
read (10,*)
do k=1,Nham
read(10,*) (eps(j,j,k), j=1,LJgroups)
end do
read (10,*)
read (10,*)
! number of different kappa's and lambda's.
num_ij= 0.5*LJgroups*(LJgroups-1)
Allocate(kappa(num_ij,Nham)) ! allocating kappa matrix.
Allocate(lambda(num_ij,Nham)) ! allocating lambda matrix.
do k=1,Nham
read (10,*) (kappa(i,k), i=1,Num_ij)
end do
read (10,*)
read (10,*)
do k=1,Nham
read (10,*) (lambda(j,k), j=1,Num_ij)
end do
read (10,*)
read (10,*) maxR
read (10,*)
read (10,*) Norm_inc
read (10,*)
read (10,*) small_inc
read (10,*)
read (10,*) change_R
close(10)
!******************* Begin Calculation ************************
! Finding center of mass of group c and let it = to (0,0,0)
call cm(Nc, Xc, Yc, Zc, Mc, XCcm, YCcm, ZCcm)
! Finding center of mass of group p
call cm(Np, Xp, Yp, Zp, Mp, XPcm, YPcm, ZPcm)
! Transform coordinates relative to center of mass of groupc.
do i=1,Nc
Xc(i)=Xc(i)-XCcm
Yc(i)=Yc(i)-YCcm
Zc(i)=Zc(i)-ZCcm
end do
do j=1,Np
Xp(j)=Xp(j)-XCcm
Yp(j)=Yp(j)-YCcm
Zp(j)=Zp(j)-ZCcm
end do
XPcm=XPcm-XCcm
YPcm=YPcm-YCcm
ZPcm=ZPcm-ZCcm
XCcm = 0
YCcm = 0
ZCcm = 0
! Correction for Epsilon and Sigma.
do h=1,Nham
do i=1,(LJgroups-1)
do j=(i+1),LJgroups
eps(i,j,h)=sqrt(eps(i,i,h)*eps(j,j,h))*(1-kappa(kappcount,h))
eps(j,i,h)=eps(i,j,h)
sig(i,j,h)=0.5*(sig(i,i,j)+sig(j,j,h))*(1-lambda(kappcount,h))
sig(j,i,h)=sig(i,j,h)
kappcount = kappcount+1
end do
end do
end do
open(20, file='rvsf.dat') ! Create output file.
r=sqrt(XPcm*XPcm+YPcm*YPcm+ZPcm*Zpcm) ! distance between groups.
increment = Small_inc ! Setting initial increment.
dt = increment
! Fix group C and move group P outwards starting from a distance='increment'
! ****************** Main loop ********************************************
do while (dt<=maxR)
if ((r>change_R).and.(num_points>1)) then
increment = norm_inc
end if
! Take steps forward at the given increment.
XP = XP - XPcm*(1-dt/r) ! calculates x coordinates of beads
YP = YP - YPcm*(1-dt/r) ! " " y " "
ZP = ZP - ZPcm*(1-dt/r) ! " " z " "
XPcm = (XPcm/r)*dt ! moves the CM of group p
YPcm = (YPcm/r)*dt
ZPcm = (ZPcm/r)*dt
dt = dt + increment
r=sqrt(XPcm*XPcm+YPcm*YPcm+ZPcm*Zpcm) ! Distance between groups.
totalF = 0.0 ! initialize the sum of f factors.
if (r<=disrot2) then
numtrials = numrot1
else if (r<=disrot3) then
numtrials = numrot2
else if (r<=disrot4) then
numtrials = numrot3
else
numtrials = numrot4
end if
! The following loop rotates each groups randomly at a fixed distance
! and calculate the energy.
do i=1,numtrials
!^^^^^^^^^^^^^^ FOR GROUP C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! Randomly choose axis to rotate:
! let 1=x axis, 2=y axis, 3=z axis
axisC =int(3*ran2(seed1)+1)
thetaC = 2*pi*ran2(seed2)
call rotate(Nc,axisC,thetaC,XCcm,YCcm,ZCcm,Xc,Yc,Zc,XcRot,YcRot,ZcRot)
!^^^^^^^^^^^^^^FOR GROUP P ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
axisP = int(3*ran2(seed3)+1)
thetaP = 2*pi*ran2(seed4)
call rotate(Np,axisP,thetaP,XPcm,YPcm,ZPcm,Xp,Yp,Zp,XpRot,YpRot,ZpRot)
!~~~~~~~~~~~ calculate energy between groups ~~~~~~~~~~~~~~~~~~~~~~~~
call ljmolecule(Nc,XcRot,YcRot,ZcRot,TYPEc,Np,XpRot,YpRot,ZpRot, &
TYPEp,Nham,LJgroups,EPS,SIG,ENERGY)
! Re-assigned the coordinates.
Xc=XcRot
Yc=YcRot
Zc=ZcRot
Xp=XpRot
Yp=YpRot
Zp=ZpRot
do k=1,Nham
Ffactor(k) = exp((-Energy(k))/(T))-1
end do
totalF=totalF+Ffactor
end do ! Num_Trials
AveFfactor = TotalF/Numtrials
write(20,*) r, AveFfactor*r*r
num_points=num_points+1
end do ! While
! ************** End Main Loop ***************************************
close(20)
LastR = dt-increment
call LocateR(fileinput,num_points,norm_inc,small_inc,MidR)
call integrate(fileinput,num_points,increment,small_inc,MidR,LastR,IntVal)
b2val = -(3*Intval)/(sig(1,1,1)*sig(1,1,1)*sig(1,1,1))
b2final = 2.094395102*0.6022*sig(1,1,1)*sig(1,1,1)*sig(1,1,1)*b2val
print*, 'The value of the second virial coefficient is : ', b2final
call result(IntVal,b2val,MidR,LastR,num_points,T,disrot1,disrot2,disrot3, &
disrot4,numrot1,numrot2,numrot3,numrot4,small_inc,norm_inc, &
Name,LJgroups,Nham,sig,eps,b2final,fileresult)
end program b2
!**************** Subroutines and functions **************
!--------------------------------------------------
subroutine cm(N, X, Y, Z, M, Xcm, Ycm, Zcm)
!--------------------------------------------------
! This subroutine calculates the center of mass of a given group.
! Inputs: N = number of beads in the group.
! X,Y,Z = coordinates of the beads.
! M = mass of the beads.
! Xcm, Ycm, Zcm = coordinates of the center of mass.
implicit none
integer, intent(in) ::N
real, dimension(N), intent(in) ::x, y, z
real, dimension(N), intent(in) ::M
real, intent(out) ::Xcm, Ycm, Zcm
integer i,j
real XM, Ym, Zm ! mass * coordinates.
real SumXM, SumYM, SumZM ! Sum of mass * coord.
real totalmass ! total mass of group.
SumXM = 0; SumYM = 0; SumZM = 0 ! Initializing.
totalmass = 0
do i=1,N
XM=x(i)*m(i)
SumXM = SumXM + XM
YM=y(i)*m(i)
SumYM = SumYM + YM
ZM=z(i)*m(i)
SumZM = SumZM + ZM
totalmass = totalmass + m(i)
end do
Xcm = SumXM/totalmass
Ycm = SumYM/totalmass
Zcm = SumZM/totalmass
return
End subroutine cm
!----------------------------------------------------------
FUNCTION ran2(idum)
!----------------------------------------------------------
INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
REAL ran2,AM,EPS,RNMX
PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, &
IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, &
NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER idum2,j,k,iv(NTAB),iy
SAVE iv,iy,idum2
DATA idum2/123456789/, iv/NTAB*0/, iy/0/
if (idum.le.0) then
idum=max(-idum,1)
idum2=idum
do j=NTAB+8,1,-1
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
if (j.le.NTAB) iv(j)=idum
enddo
iy=iv(1)
endif
k=idum/IQ1
idum=IA1*(idum-k*IQ1)-k*IR1
if (idum.lt.0) idum=idum+IM1
k=idum2/IQ2
idum2=IA2*(idum2-k*IQ2)-k*IR2
if (idum2.lt.0) idum2=idum2+IM2
j=1+iy/NDIV
iy=iv(j)-idum2
iv(j)=idum
if(iy.lt.1)iy=iy+IMM1
ran2=min(AM*iy,RNMX)
return
END
!------------------------------------------------------------------------
subroutine rotate(N,axis,theta,Xcm,Ycm,Zcm,Xold,Yold,Zold,Xnew,Ynew,Znew)
!------------------------------------------------------------------------
! This subroutine rotates the molecule for a given theta around a
! specified axis.
implicit none
integer, intent(in) :: N ! number of beads in a group.
integer, intent(in) :: axis
double precision :: theta
real, dimension(N), intent(in) :: Xold,Yold,Zold
real, dimension(N), intent(out) :: Xnew,Ynew,Znew
real, dimension(N) :: Xtemp,Ytemp,Ztemp
real, intent(in) :: Xcm,Ycm,Zcm
integer k
Xtemp = Xold-Xcm
Ytemp = Yold-Ycm
Ztemp = Zold-Zcm
if (axis==1) then
do k=1,N
Ynew(k)=Ytemp(k)*cos(theta) - Ztemp(k)*sin(theta)+Ycm
Znew(k)=Ytemp(k)*sin(theta) + Ztemp(k)*cos(theta)+Zcm
Xnew(k)=Xtemp(k)+Xcm
enddo
else if (axis==2) then
do k=1,N
Xnew(k)=Xtemp(k)*cos(theta) - Ztemp(k)*sin(theta)+Xcm
Znew(k)=Xtemp(k)*sin(theta) + Ztemp(k)*cos(theta)+Zcm
Ynew(k)=Ytemp(k)+Ycm
enddo
else
do k=1,N
Xnew(k)=Xtemp(k)*cos(theta) - Ytemp(k)*sin(theta)+Xcm
Ynew(k)=Xtemp(k)*sin(theta) + Ytemp(k)*cos(theta)+Ycm
Znew(k)=Ztemp(k)+Zcm
enddo
end if
return
end subroutine rotate
!----------------------------------------------------------------
subroutine ljmolecule(Nc,Xc,Yc,Zc,TYPEc,Np,Xp,Yp,Zp, &
TYPEp,Nham,Nljgrs,EPS,SIG,ENERGY)
!----------------------------------------------------------------
! This subroutine calculates the energy between two LJ groups in space.
! Nc,Np = number of beads in groups c and p.
! Xc,Yc,Zc,Xp,Yp,Zp = coordinate of the beads
! Nham = number of hamiltonians.
! Nljgrs = number of LJ groups.
! EPS,SIG = Rank 3 matrix for epsilon and sigma.
! Energy = energy of interaction between groups for each
! hamiltonian.
integer, intent(in) :: Nc,Np
integer, intent(in) :: Nham
integer, dimension(Nc), intent(in) :: TYPEc
integer, dimension(Np), intent(in) :: TYPEp
integer, intent(in) :: Nljgrs
real, dimension(Nc), intent(in) :: Xc, Yc, Zc
real, dimension(Np), intent(in) :: Xp, Yp, Zp
real, dimension(Nljgrs,Nljgrs,Nham) :: EPS, SIG
real, dimension(Nham), intent(out) :: Energy
integer :: i, j, h
real :: xij, yij, zij
real :: sigor2, sigor6
Energy = 0.0
do i=1,Nc
do j=1,Np
xij = abs(Xp(j) - Xc(i))
yij = abs(Yp(j) - Yc(i))
zij = abs(Zp(j) - Zc(i))
rij2 = xij*xij + yij*yij + zij*zij
do h = 1,Nham
sigor2=(SIG(TYPEp(j),TYPEc(i),h) * &
SIG(TYPEp(j),TYPEc(i),h))/rij2
sigor6=sigor2*sigor2*sigor2
ENERGY(h) =Energy(h)+4*(sigor6*sigor6-sigor6)* &
EPS(TYPEp(j),TYPEc(i),h)
end do
end do
end do
return
end subroutine ljmolecule
!----------------------------------------------------------------------------
subroutine LocateR(Inputfile,num_points,Norm_Inc,Min_Inc, &
MidR)
!----------------------------------------------------------------------------
! This subroutine sorts out the unwanted data points and
! locates the position where the increments decrease or increase
! and the number of segments that have the smaller increment.
implicit none
character*20, intent(in) :: InputFile ! Main input file
integer p
logical first
real, intent(out) :: MidR
real, intent(in) :: Norm_Inc
real, intent(in) :: Min_Inc
real r1,r2,r3
integer num_points
real diff1,diff2
real avef2, avef3
open(20,file=Inputfile)
r1 = 0.0
first = .true.
read(20,*) r2, avef2
Do p=2,Num_Points
if (first) then
read(20,*) r3, avef3
diff1 = r2-r1
diff2 = r3-r2
if (abs(diff2-diff1)>0.0001) then
MidR = r2
first= .not. first
end if
r1 = r2
r2 = r3
avef2 = avef3
end if
End do
close(20)
return
end subroutine LocateR
!---------------------------------------------------------------
subroutine integrate(Inputfile,Num_points,Norm_inc,Small_inc, &
MidR,LastR,IntVal)
!---------------------------------------------------------------
! This subroutine integrates f(r) using data from the output
! file called rvsf.dat. The integration methods used are Simpson's
! 1/3 and 3/8 where applicable.
implicit none
character*20, intent(in) :: InputFile ! Main input file
real, intent(out) :: IntVal ! value of integral.
real, intent(in) :: MidR, LastR
real, intent(in) :: Norm_inc ! value of increment.
real, intent(in) :: Small_inc
real, dimension(Num_points+1) :: f
real, dimension(Num_points+1) :: r
real h, hfuture ! Length of a segment.
integer n ! number of segments.
integer num_points
integer i,j,k ! counters.
!***************** Read Data from Output file *****************
open(40,file=Inputfile)
! Find number of segments.
n = nint(MidR/small_inc)+nint((LastR-MidR)/Norm_inc)
IntVal = 0.0
r(0) = 0.0
f(0) = 0.0
do i=1,n
read(40,*) r(i), f(i)
end do
h = r(1)-r(0)
k = 1
do j=1,n
hfuture = r(j+1) - r(j)
if (abs(h-hfuture)<0.0001) then
if (k==3) then
IntVal = IntVal + 2.0*h*(f(j-1)+4.0*f(j-2)+f(j-3))/6.0
k = k-1
else
k = k+1
end if
else if (k==1) then
IntVal = IntVal + h*(f(j)+f(j-1))/2.0
else if (k==2) then
IntVal = IntVal +2.0*h*(f(j)+4.0*f(j-1)+f(j-2))/6.0
k = 1
else
IntVal = IntVal + 3.0*h*(f(j)+3.0*(f(j-1)+f(j-2))+f(j-3))/8.0
k = 1
end if
h = hfuture
end do
close(40)
return
end subroutine integrate
! ------------------------------------------------------------------------------
subroutine result(IntVal,b2val,MidR,LastR,num_points,T,disrot1,disrot2, &
disrot3,disrot4,numrot1,numrot2,numrot3,numrot4, &
small_inc,norm_inc,Name,LJgroups,Nham,sig,eps,b2final,fileresult)
! ------------------------------------------------------------------------------
real, intent(in) :: Intval
real, intent(in) :: b2val
real, intent(in) :: MidR,LastR
real, intent(in) :: T
real, intent(in) :: disrot1,disrot2,disrot3,disrot4
integer, intent(in) :: Nham, LJgroups
real, dimension(LJgroups,LJgroups,Nham), intent(in) :: sig,eps
real, intent(in) :: norm_inc, small_inc
integer, intent(in) :: numrot1,numrot2,numrot3,numrot4
character*20, intent(out) :: fileresult
character*20, intent(in) :: Name
real, intent(in) :: b2final
open(6,file=fileresult)
write(6,*) 'Results file for second virial coefficient calculation'
write(6,*)
write(6,*) 'Calculation for molecule type : ', Name
write(6,*)
write(6,*) 'Temperature = ',T, 'K'
write(6,*)
write(6,*) 'Epsilon/k (K) = ', eps
write(6,*)
write(6,*) 'Sigma (Angstrom) = ', sig
write(6,*)
write(6,*) 'Maximum distance evaluated is at r = ', LastR, ' Angstrom '
write(6,*)
write(6,*) 'Number of rotations from ',disrot1, 'to',disrot2,' = ',numrot1
write(6,*)
write(6,*) 'Number of rotations from ',disrot2, 'to',disrot3,' = ',numrot2
write(6,*)
write(6,*) 'Number of rotations from ',disrot3, 'to',disrot4,' = ',numrot3
write(6,*)
write(6,*) 'Number of rotations from ',disrot4, 'to',LastR,' = ',numrot4
write(6,*)
write(6,*) 'Number of data points integrated: ', num_points
write(6,*)
write(6,*) 'Small increment = ', small_inc, ' Angstrom starts at 0 and changes'
write(6,*) 'to increment = ', norm_inc , 'at r = ', MidR, ' Angstrom'
write(6,*)
write(6,*) 'The value of the integral is ', IntVal
write(6,*)
write(6,*) 'The reduced value (B*) of the second virial coefficient is ', b2val
write(6,*)
write(6,*) 'The value of the second virial coefficient is ', b2final, ' cm^3/mol'
close(6)
return
end subroutine result
Back to the 