// First Displacement of particle and Acceptance of Move
subroutine displace(success, ihs)
Use sp
Use xyz
Use randmod
Use moves
Use tmmod
implicit none
logical :: success
integer :: ihs,i
integer :: mol
real,dimension(1:5) :: xold, yold, zold
real,dimension(1:5) :: xnew, ynew, znew
real :: dx, dy, dz, xo,yo,zo,xn,yn,zn
real :: emffnew, emffold, ioneffold, ioneffnew, ioneffold1, ioneffnew1
real :: emwfnew, emwfold, ionewfnew, ionewfold, ionewfnew1, ionewfold1
real :: de, dlnpsi
success = .false.
emffnew = 0
emffold = 0
emwfnew = 0
emwfold = 0
ioneffold = 0
ioneffold1 = 0
ioneffnew = 0
ioneffnew1 = 0
ionewfold = 0
ionewfold1 = 0
ionewfnew = 0
ionewfnew1 = 0
pacc = 0.0
if( Nmol == 0 ) return
ihs = 1
if( ran2(Seed) > 0.33 ) ihs = 2
if( ihs == 1 ) then
dx = 0.0
dy = 0.0
dz = ( 2.0 * ran2(Seed) - 1.0 ) * dh * bh
else
dx = ( 2.0 * ran2(Seed) - 1.0 ) * ds * bs
dy = ( 2.0 * ran2(Seed) - 1.0 ) * ds * bs
dz = 0.0
end if
mol = int( Nmol * ran2(Seed) ) + 1
if (mol .eq. Nmol) then
xo = Xr
yo = yr
zo = zr
xn = xo + dx
yn = yo + dy
zn = zo + dz
if( xn > bs ) xn = xn - bs * aint( xn / bs )
if( yn > bs ) yn = yn - bs * aint( yn / bs )
if( nwalls == 0 .and. zn > bh ) zn = zn - bh * aint( zn / bh )
if( nwalls > 0 .and. zn > bh ) return
if( xn < 0.0 ) xn = xn - bs * aint( xn / bs - 1.0 )
if( yn < 0.0 ) yn = yn - bs * aint( yn / bs - 1.0 )
if( nwalls == 0 .and. zn < 0.0 ) zn = zn - bh * aint( zn / bh - 1.0 )
if( nwalls > 0 .and. zn < 0.0 ) return
call ljmoleculeion(ioneffold,ionewfold)
call ljmoleculeion1(ioneffold1,ionewfold1)
ioneffold = (1-lambda) * ioneffold + (lambda * ioneffold1)
Xr = xn
Yr = yn
Zr = zn
call ljmoleculeion(ioneffnew,ionewfnew)
call ljmoleculeion1(ioneffnew1,ionewfnew1)
ioneffnew = (1-lambda) * ioneffnew + (lambda * ioneffnew1)
de = (ioneffnew + ionewfnew ) - (ioneffold + ionewfold)
else
xold(:) = X(mol,:)
yold(:) = Y(mol,:)
zold(:) = Z(mol,:)
xnew = xold + dx
ynew = yold + dy
znew = zold + dz
do i=1,5
if( xnew(i) > bs ) xnew(i) = xnew(i) - bs * aint( xnew(i) / bs )
if( ynew(i) > bs ) ynew(i) = ynew(i) - bs * aint( ynew(i) / bs )
if( nwalls == 0 .and. znew(i) > bh ) znew(i) = znew(i) - bh * aint( znew(i) / bh )
if( nwalls > 0 .and. znew(i) > bh ) return
if( xnew(i) < 0.0 ) xnew(i) = xnew(i) - bs * aint( xnew(i) / bs - 1.0 )
if( ynew(i) < 0.0 ) ynew(i) = ynew(i) - bs * aint( ynew(i) / bs - 1.0 )
if( nwalls == 0 .and. znew(i) < 0.0 ) znew(i) = znew(i) - bh * aint( znew(i) / bh - 1.0 )
if( nwalls > 0 .and. znew(i) < 0.0 ) return
end do
call ljmolecule(mol, emffold, emwfold)
X(mol,:) = xnew(:)
Y(mol,:) = ynew(:)
Z(mol,:) = znew(:)
call ljmolecule(mol, emffnew, emwfnew)
de = ( emffnew + emwfnew ) - ( emffold + emwfold )
endif
dlnpsi = - beta * de
pacc = -beta * de
if ( de > 2160.0) then
!write(*,*)' problem is in displaces stop stop stop stopstop*************************',mol
write(*,*)' change',ihs,dx,dy,dz,de
write(*,*) ' energies', emffnew,emffold
! write(*,*)' xold', xold
! write(*,*)' yold', yold
!write(*,*)' zold', zold
!write(*,*)' xnew', xnew
!write(*,*)' ynew', ynew
!write(*,*)' znew', znew
end if
if( pacc < 0.0 ) then
pacc = exp(pacc)
else
pacc = 1.0
end if
if( log( ran2(Seed) ) < dlnpsi ) then
success = .True.
!ebff = ebff + emffnew - emffold + ioneffnew - ioneffold
ebff=ebff+de
ebwf = ebwf + emwfnew - emwfold + ionewfnew - ionewfold
else
if (mol .eq. Nmol) then
xr = xo
yr = yo
zr = zo
else
X(mol,:) = xold(:)
Y(mol,:) = yold(:)
Z(mol,:) = zold(:)
endif
end if
return
end subroutine displace
//LJ Potential calculation
subroutine ljinteract(eff,ewf)
Use xyz
Use sp
Use lj
implicit none
real :: eff,uc1, uc2, uc3, ewf
integer :: i, j
real :: xij, yij, zij
real :: dx,dy,dz
real :: xi, yi, zi
real :: ioneff0, ionewf0, ioneff1, ionewf1, ioneff, ionewf
real :: rij2, rj2, xcmij, ycmij, zcmij, roij2
real, external :: uff, uwf, ucharge, usflj
!****************************************************************
eff = 0.0
ewf = 0.0
if( Nmol == 0 ) return
do i = 1, Nmol - 2
do j = i+1, Nmol-1
xi = X(i,5)
yi = Y(i,5)
zi = Z(i,5)
xij = abs( X(j,5) - xi )
yij = abs( Y(j,5) - yi )
zij = abs( Z(j,5) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rij2 = xij*xij + yij*yij + zij*zij
xi = X(i,1)
yi = Y(i,1)
zi = Z(i,1)
xij = abs( X(j,1) - xi )
yij = abs( Y(j,1) - yi )
zij = abs( Z(j,1) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
roij2 = xij*xij + yij*yij + zij*zij ! Distance between O-O atoms
if( roij2 < rcutsq ) then
eff = eff + usflj(roij2)
end if
!m1
xi = X(i,2)
yi = Y(i,2)
zi = Z(i,2)
!m1 m2
xij = abs( X(j,2) - xi )
yij = abs( Y(j,2) - yi )
zij = abs( Z(j,2) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,1,1)
!m1 h21
xij = abs( X(j,3) - xi )
yij = abs( Y(j,3) - yi )
zij = abs( Z(j,3) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,1,2)
!m1 h22
xij = abs( X(j,4) - xi )
yij = abs( Y(j,4) - yi )
zij = abs( Z(j,4) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,1,2)
!H11
xi = X(i,3)
yi = Y(i,3)
zi = Z(i,3)
!H11 m2
xij = abs( X(j,2) - xi )
yij = abs( Y(j,2) - yi )
zij = abs( Z(j,2) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,1)
!H11 H21
xij = abs( X(j,3) - xi )
yij = abs( Y(j,3) - yi )
zij = abs( Z(j,3) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
!H11 H22
xij = abs( X(j,4) - xi )
yij = abs( Y(j,4) - yi )
zij = abs( Z(j,4) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
!H12
xi = X(i,4)
yi = Y(i,4)
zi = Z(i,4)
!H12, M2
xij = abs( X(j,2) - xi )
yij = abs( Y(j,2) - yi )
zij = abs( Z(j,2) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,1)
!H12 H21
xij = abs( X(j,3) - xi )
yij = abs( Y(j,3) - yi )
zij = abs( Z(j,3) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
!H12 H22
xij = abs( X(j,4) - xi )
yij = abs( Y(j,4) - yi )
zij = abs( Z(j,4) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
end do
if( nwalls > 0 ) ewf = ewf + uwf(zi)
if( nwalls == 2) ewf = ewf + uwf(bh - zi)
end do
zi = Z(Nmol,5)
if( nwalls > 0 ) ewf = ewf + uwf(zi)
if( nwalls == 2) ewf = ewf + uwf(bh - zi)
! ** To calculate energy of water + Metal ion
call ljmoleculeion(ioneff0,ionewf0)
call ljmoleculeion1(ioneff1,ionewf1)
if(mutation) then
ioneff=(1-lambda)*ioneff0 + (lambda*ioneff1)
ionewf=(1-lambda)*ionewf0 + (lambda*ionewf1)
else
ioneff = ioneff0
ionewf = ionewf0
end if
eff = eff + ioneff
ewf = ewf + ionewf
return
end subroutine ljinteract
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ljmoleculeion(ioneff,ionewf)
Use xyz
Use sp
Use lj
implicit none
real :: ioneff, ioneff1, ionewf,uc1, uc2, uc3
integer :: i, j, etotal
real :: xij, yij, zij
real :: dx, dy, dz
real :: xi, yi, zi
real :: rij2, rj2, xcmij, ycmij, zcmij, roij2
real, external :: uff, uwf, ucharge, uljion, uljion1,uchargeion
!****************************************************************
ioneff = 0.0
ionewf = 0.0
if( Nmol == 0 ) return
do j = 1, Nmol-1
xij = abs( X(j,5) - Xr )
yij = abs( Y(j,5) - Yr )
zij = abs( Z(j,5) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rij2 = xij*xij + yij*yij + zij*zij
!ion-O vander wall interaction
xij = abs( X(j,1) - Xr )
yij = abs( Y(j,1) - Yr )
zij = abs( Z(j,1) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
roij2 = xij*xij + yij*yij + zij*zij
if( roij2 < rcutsq ) then
ioneff = ioneff + uljion(roij2)
end if
!ion-M coulombic interaction
xij = abs( X(j,2) - Xr )
yij = abs( Y(j,2) - Yr )
zij = abs( Z(j,2) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
ioneff = ioneff + ucharge(rj2,roij2,3,1)
!ion-H1 coulombic interaction
xij = abs( X(j,3) - Xr )
yij = abs( Y(j,3) - Yr )
zij = abs( Z(j,3) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
ioneff = ioneff + ucharge(rj2,roij2,3,2)
!ion-H2 coulombic interaction
xij = abs( X(j,4) - Xr )
yij = abs( Y(j,4) - Yr )
zij = abs( Z(j,4) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
ioneff = ioneff + ucharge(rj2,roij2,3,2)
end do
if( nwalls > 0 ) ionewf = ionewf + uwf(Zr)
if( nwalls == 2) ionewf = ionewf + uwf(bh - Zr)
return
end subroutine ljmoleculeion
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ljmoleculeion1(ioneff1,ionewf1)
Use xyz
Use sp
Use lj
implicit none
real :: ioneff1, ionewf1, uc1, uc2, uc3
integer :: i, j, etotal
real :: xij, yij, zij
real :: dx, dy, dz
real :: xi, yi, zi
real :: rij2, rj2, xcmij, ycmij, zcmij, roij2
real, external :: uff, uwf, ucharge, uljion1, uchargeion
!****************************************************************
ioneff1 = 0.0
ionewf1 = 0.0
if( Nmol == 0 ) return
do j = 1, Nmol-1
xij = abs( X(j,5) - Xr )
yij = abs( Y(j,5) - Yr )
zij = abs( Z(j,5) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rij2 = xij*xij + yij*yij + zij*zij
!ion-O vander wall interaction
xij = abs( X(j,1) - Xr )
yij = abs( Y(j,1) - Yr )
zij = abs( Z(j,1) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
roij2 = xij*xij + yij*yij + zij*zij
if( roij2 < rcutsq ) then
ioneff1 = ioneff1 + uljion1(roij2)
end if
!ion-M coulombic interaction
xij = abs( X(j,2) - Xr )
yij = abs( Y(j,2) - Yr )
zij = abs( Z(j,2) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
ioneff1 = ioneff1 + ucharge(rj2,roij2,3,1)
!ion-H1 coulombic interaction
xij = abs( X(j,3) - Xr )
yij = abs( Y(j,3) - Yr )
zij = abs( Z(j,3) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
ioneff1 = ioneff1 + ucharge(rj2,roij2,3,2)
!ion-H2 coulombic interaction
xij = abs( X(j,4) - Xr )
yij = abs( Y(j,4) - Yr )
zij = abs( Z(j,4) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
ioneff1 = ioneff1 + ucharge(rj2,roij2,3,2)
end do
if( nwalls > 0 ) ionewf1 = ionewf1 + uwf(Zr)
if( nwalls == 2) ionewf1 = ionewf1 + uwf(bh - Zr)
return
end subroutine ljmoleculeion1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ljmolecule(mol, eff, ewf)
Use xyz
Use sp
Use lj
implicit none
real :: eff, uc1, uc2, uc3, ewf
real ::eff0, eff1
integer :: i, j, mol
real :: xij, yij, zij
real :: dx, dy, dz
!real :: xt, yt, zt
real :: xi, yi, zi
real :: rij2, rj2, xcmij, ycmij, zcmij, roij2
real, external :: uff, uwf, ucharge, usflj, uljion, uljion1, uchargeion
!****************************************************************
eff = 0.0
ewf = 0.0
eff0= 0.0
eff1= 0.0
do j = 1, Nmol-1
if( j == mol ) cycle
xi = X(mol,5)
yi = Y(mol,5)
zi = Z(mol,5)
xij = abs( X(j,5) - xi )
yij = abs( Y(j,5) - yi )
zij = abs( Z(j,5) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rij2 = xij*xij + yij*yij + zij*zij
xi = X(mol,1)
yi = Y(mol,1)
zi = Z(mol,1)
xij = abs( X(j,1) - xi )
yij = abs( Y(j,1) - yi )
zij = abs( Z(j,1) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
roij2 = xij*xij + yij*yij + zij*zij
if( roij2 < rcutsq ) then
eff = eff + usflj(roij2)
end if
!m1
xi = X(mol,2)
yi = Y(mol,2)
zi = Z(mol,2)
!m1 m2
xij = abs( X(j,2) - xi )
yij = abs( Y(j,2) - yi )
zij = abs( Z(j,2) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,1,1)
!m1 H21
xij = abs( X(j,3) - xi )
yij = abs( Y(j,3) - yi )
zij = abs( Z(j,3) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,1,2)
!m1 H22
xij = abs( X(j,4) - xi )
yij = abs( Y(j,4) - yi )
zij = abs( Z(j,4) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,1,2)
!H11
xi = X(mol,3)
yi = Y(mol,3)
zi = Z(mol,3)
!H11 m2
xij = abs( X(j,2) - xi )
yij = abs( Y(j,2) - yi )
zij = abs( Z(j,2) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,1)
!H11 H21
xij = abs( X(j,3) - xi )
yij = abs( Y(j,3) - yi )
zij = abs( Z(j,3) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
!H11 H22
xij = abs( X(j,4) - xi )
yij = abs( Y(j,4) - yi )
zij = abs( Z(j,4) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
!H12
xi = X(mol,4)
yi = Y(mol,4)
zi = Z(mol,4)
!H12, M2
xij = abs( X(j,2) - xi )
yij = abs( Y(j,2) - yi )
zij = abs( Z(j,2) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,1)
!H12 H21
xij = abs( X(j,3) - xi )
yij = abs( Y(j,3) - yi )
zij = abs( Z(j,3) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
!H12 H22
xij = abs( X(j,4) - xi )
yij = abs( Y(j,4) - yi )
zij = abs( Z(j,4) - zi )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff = eff + ucharge(rj2,roij2,2,2)
end do
! Code for Interaction of water with Ion(r)
xij = abs( X(mol,5) - Xr )
yij = abs( Y(mol,5) - Yr )
zij = abs( Z(mol,5) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rij2 = xij*xij + yij*yij + zij*zij
! Ion-O vander wall interaction
xij = abs( X(mol,1) - Xr )
yij = abs( Y(mol,1) - Yr )
zij = abs( Z(mol,1) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
roij2 = xij*xij + yij*yij + zij*zij
if( roij2 < rcutsq ) then
eff0 = uljion(roij2)
eff1 = uljion1(roij2)
end if
! Ion-M coulombic interaction
xij = abs( X(mol,2) - Xr )
yij = abs( Y(mol,2) - Yr )
zij = abs( Z(mol,2) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff0 = eff0 + ucharge(rj2,roij2,3,1)
eff1 = eff1 + ucharge(rj2,roij2,3,1)
!ion-H1 coulombic interaction
xij = abs( X(mol,3) - Xr )
yij = abs( Y(mol,3) - Yr )
zij = abs( Z(mol,3) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff0 = eff0 + ucharge(rj2,roij2,3,2)
eff1 = eff1 + ucharge(rj2,roij2,3,2)
!ion-H2 coulombic interaction
xij = abs( X(mol,4) - Xr )
yij = abs( Y(mol,4) - Yr )
zij = abs( Z(mol,4) - Zr )
if( xij > bs - xij ) xij = xij - bs
if( yij > bs - yij ) yij = yij - bs
if( zij > bh - zij ) zij = zij - bh
rj2 = xij*xij + yij*yij + zij*zij
eff0 = eff0 + ucharge(rj2,roij2,3,2)
eff1 = eff1 + ucharge(rj2,roij2,3,2)
if( nwalls > 0 ) ewf = ewf + uwf(zi)
if( nwalls == 2) ewf = ewf + uwf(bh - zi)
eff = eff + (1-lambda)*eff0+(lambda*eff1)
return
end subroutine ljmolecule
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function usflj(rij2)
Use lj
implicit none
real :: usflj, rij2, rij
real, external :: ulj
select case(iljpot)
case(1)
usflj = ulj(rij2)
case(2)
usflj = ulj(rij2) - uljcut
case(3)
rij = sqrt(rij2)
usflj = ulj(rij2) - uljcut - (rij - rcut) * upljcut
end select
return
end function usflj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function ulj(rij2)
Use lj
implicit none
real :: ulj, rij2, r2, r6, r12
r2 = 1.0 / rij2
r6 = r2 * r2 * r2
r12 = r6 * r6
ulj = 4.0 * (r12 - r6)
return
end function ulj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function uljion(rij2)
Use lj
implicit none
real :: uljion, rij2, r2, r6, r12
real :: sigmar2, sigmar6, sigmar12
sigmar2 = ionsigmar * ionsigmar
sigmar6 = sigmar2 * sigmar2 * sigmar2
sigmar12= sigmar6 * sigmar6
r2 = 1.0 / rij2
r6 = r2 * r2 * r2
r12 = r6 * r6
uljion = 4.0 *ionepsr* (r12*sigmar12 - r6*sigmar6)
return
end function uljion
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function uljion1(rij2)
Use lj
implicit none
real :: uljion1, rij2, r2, r6, r12
real :: sigmam2, sigmam6, sigmam12
sigmam2 = ionsigmam * ionsigmam
sigmam6 = sigmam2 * sigmam2 * sigmam2
sigmam12= sigmam6 * sigmam6
r2 = 1.0 / rij2
r6 = r2 * r2 * r2
r12 = r6 * r6
uljion1 = 4.0 *ionepsm* (r12*sigmam12 - r6*sigmam6)
return
end function uljion1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function uplj(rij)
Use lj
implicit none
real :: uplj, rij, r2, r6, r12
r2 = 1.0 / (rij * rij)
r6 = r2 * r2 * r2
r12 = r6 * r6
uplj = -24.0 * (2.0*r12 - r6) / rij
return
end function uplj
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!**********************
function ucharge(r2,rm2,i,j)
!Use tip4p
Use lj
implicit none
real :: ucharge,r2,rm2
integer :: i,j,k
!if (r2 .le. 0.16) then
! ucharge=1.0E100
!else
if(rm2 <= rcutsq) then
ucharge = rb * charge(i)*charge(j)/sqrt(r2)
else
ucharge=0.0
end if
end function ucharge
!*****************************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!***************
function uchargeion(r2,rm2,i,j)
Use lj
implicit none
real :: uchargeion,r2,rm2
integer :: i,j,k
if(rm2 <= rcutsq) then
uchargeion = 215.8323 * charge(i)*charge(j)/sqrt(r2)
else
uchargeion=0.0
end if
end function uchargeion
!*****************************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function uwf(z)
Use lj
implicit none
real :: uwf, z, ar, aa, ap, r, r2, r3, r4, r9
select case(isfpot)
case(1)
r = wsig / z
r3 = r * r * r
r9 = r3 * r3 * r3
aa = r3
ar = 2.0 / 15.0 * r9
uwf = weps * (ar - aa)
case(2)
r = wsig / z
r2 = r * r
r4 = r2 * r2
aa = r4
ar = 2.0 / 5.0 * r4 * r4 * r2
ap = z + pot2
ap = pot1 / (ap * ap * ap)
uwf = weps * (ar - aa - ap)
end select
return
end function uwf
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
No comments:
Post a Comment