subroutine Volume_NPT( Nmol, Nlj, Nion, MaxMol, Xlj, Ylj, Zlj, TYPElj, & Xion, Yion, Zion, TYPEion, BoxSize, Constant, dvol, & Pressure, LENGTHlj, LENGTHion, STARTlj, STARTion, & NTemp, BETA, Nham, LNWSTATE, XLRCORR, ELRCORR,& Nljgrs, EPS, SIG, CP, ALP, RMAX, MASSlj, Niongrs, & CHARGE, MASSion, Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, & EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, & ULJBOX, ULJLR, UREAL, UFOURIER, USURF, USELF_CH, & USELF_MOL, Increase, Success, TEST_PRESS, Seed ) implicit none ! This subroutine is used to attempt a volume change of a simulation box. ! Nmol is the number of molecules in the simulation box. ! Nlj is the number of LJ beads in the simulation box. ! Nion is the number of ionic beads in the simulation box. integer, intent(in) :: Nmol integer, intent(in) :: Nlj integer, intent(in) :: Nion ! MaxMol is the maximum number of molecules. integer :: MaxMol ! Xlj, Ylj, and Zlj are the coordinates of the LJ beads. ! Xion, Yion, and Zion are the coordinates of the ionic beads. real, dimension( Nlj ), intent(inout) :: Xlj, Ylj, Zlj real, dimension( Nion ), intent(inout) :: Xion, Yion, Zion ! TYPElj contains the group identity of each LJ bead. ! TYPEion contains the group identity of each ionic bead. integer, dimension( Nlj ), intent(in) :: TYPElj integer, dimension( Nion ), intent(in) :: TYPEion ! BoxSize is the length of the simulation box. real, intent(inout) :: BoxSize ! Costant is a logical indicating whether the the Volume changes are ! done with a constant DV or variable DV, constant DV ==> Constant = .True. ! dvol is the constant DV or the maximum DV. logical, intent(in) :: Constant real, intent(in) :: dvol ! Pressure is the pressure for an NPT run. Set Pressure = 0 for a constant ! volume run in which this subroutine is used to calculate the pressure. ! Pressure should have units of bar. real, intent(in) :: Pressure ! LENGTHlj contains the number of LJ beads in each molecule. ! LENGTHion contains the number of ionic beads in each molecule. ! STARTlj contains the starting LJ bead number for each molecule. ! STARTion contains the starting ionic bead number for each molecule. integer, dimension( Nmol ), intent(in) :: LENGTHlj, LENGTHion integer, dimension( Nmol ), intent(in) :: STARTlj, STARTion ! NTemp is the number of temperatures being used. integer, intent(in) :: NTemp ! BETA contains the reciprical temperature. real, dimension(NTemp,1), intent(in) :: BETA ! Nham is the number of hamiltonians being used. integer, intent(in) :: Nham ! LNWSTATE contains the weight of each temperature and hamiltonian. real, dimension(NTemp, Nham), intent(inout) :: LNWSTATE ! XLRCORR and ELRCORR contain parameters for the long range LJ energy. real, dimension(605), intent(in) :: XLRCORR, ELRCORR ! Nljgrs is the number of LJ groups in the system. ! EPS is a rank 3 array containing the eps_ij parameters for each hamiltonian. ! SIG is a rank 3 array containing the sigma_ij parameters for each hamiltonian. ! CP contains the C_ij parameters for each hamiltonian. ! ALP contains the alpha_ij parameters for each hamiltonian. ! RMAX contains the Rmax_ij parameters for each hamiltonian. integer, intent(in) :: Nljgrs real, dimension(Nljgrs, Nljgrs, Nham), intent(in) :: EPS, SIG, CP, ALP, RMAX ! MASSlj contains the mass of the LJ groups. real, dimension(Nljgrs), intent(in) :: MASSlj ! Niongrs is the number of ionic groups in the system. ! CHARGE is a rank 2 array containing the charge of group i for each hamiltonian. integer, intent(in) :: Niongrs real, dimension(Niongrs, Nham), intent(in) :: CHARGE ! MASSion contains the mass of the ionic groups. real, dimension(Niongrs), intent(in) :: MASSion ! Alpha is an Ewald sum parameter, Alpha = kappa * L, for kappa in A + T. real, intent(in) :: Alpha ! Kmax is an Ewald sum parameter. ! Nkvec is the number of k-vectors used in the Fourier sum. ! KX, KY, KZ contain the vector identity of the Nkvec vectors. ! CONST contains the constant part of the Fourier summation for a given Nkvec. integer, intent(in) :: Kmax integer, intent(in) :: Nkvec integer, dimension(Nkvec), intent(in) :: KX, KY, KZ real, dimension(Nkvec), intent(in) :: CONST ! EXPX contains the value of exp( i*kx*x ) for a given kx and ion. complex, dimension(0:Kmax, Nion), intent(inout) :: EXPX complex, dimension(-Kmax:Kmax,Nion), intent(inout) :: EXPY, EXPZ ! SUMQEXPV contains the summation of qi*exp(i*(kx*x + ky*y + kz*z)) ! for a given k-vector and hamiltonian. complex, dimension(Nkvec, Nham), intent(inout) :: SUMQEXPV ! SUMQX is the summation of qi * xi for all ions in the box. real, dimension(Nham), intent(inout) :: SUMQX, SUMQY, SUMQZ ! ULJ is the LJ energy of the system without the long range correction. ! UFOURIER is the coulombic fourier energy of the system. ! UREAL is the coulombic real energy of the system. ! USURF is the coulombic surface energy of the system. real, dimension(Nham), intent(inout) :: ULJBOX, ULJLR real, dimension(Nham), intent(inout) :: UREAL, UFOURIER real, dimension(Nham), intent(inout) :: USURF, USELF_CH real, dimension(MaxMol,Nham), intent(inout) :: USELF_MOL ! Increase is a flag to indicate whether the volume was increased or decreased. ! Increase = 1 for an increase. ! Increase = 2 for a decrease. integer, intent(out) :: Increase ! Success is a logical indicating whether the move was successful or not. logical, intent(out) :: Success ! TEST_PRESS is the pressure calculated from the move. real, dimension(NTemp, Nham), intent(out) :: TEST_PRESS ! Seed is the current random number generator seed value. integer, intent(inout) :: Seed real, external :: ran2 ! Local Variables. integer :: lenlj, stlj, endlj integer :: lenion, stion, endion integer, dimension(Nljgrs) :: NGROUPS integer, dimension(Nlj + Nion) :: temp4 integer :: h, j, k logical :: Ions = .False. real :: CoulCombo real :: Largest, dV real :: PiRatio, LnPiRatio real :: xcom, ycom, zcom real :: BoxSize_new, BoxRatio real, dimension(Nlj) :: Xlj_new, Ylj_new, Zlj_new real, dimension(Nlj + Nion) :: temp1, temp2, temp3 real, dimension(Nljgrs + Niongrs) :: temp5 real, dimension(Nion) :: Xion_new, Yion_new, Zion_new real, dimension(NTemp, Nham) :: BETA_HAM real, dimension(Nham) :: ULJBOX_new, ULJLR_new real, dimension(Nham) :: UREAL_new, UFOURIER_new real, dimension(Nham) :: USURF_new, USELF_CH_new real, dimension(MaxMol,Nham) :: USELF_MOL_new real, dimension(Nham) :: SUMQX_NEW, SUMQY_NEW, SUMQZ_NEW real, dimension(Nham) :: ZERO real, dimension(Nham,1) :: dU real, dimension(Nljgrs,Nljgrs,Nham) :: EPS_TMP real, parameter :: Pi = 3.14159265359 real, parameter :: ec = 1.60217733e-19 real, parameter :: eps0 = 8.854187817e-12 real, parameter :: kB = 1.380658e-23 complex, dimension(0:Kmax , Nion) :: EXPX_NEW complex, dimension(-Kmax:Kmax , Nion) :: EXPY_NEW complex, dimension(-Kmax:Kmax , Nion) :: EXPZ_NEW complex, dimension(Nkvec, Nham) :: SUMQEXPV_NEW complex, dimension(0:Kmax , Nion) :: CZERO1 complex, dimension(-Kmax:Kmax , Nion) :: CZERO2 complex, dimension(Nkvec, Nham) :: CZERO3 Success = .False. dV = dvol if( .NOT. Constant ) then dV = dV * ran2(Seed) * BoxSize ** 3.0 end if if( ran2(Seed) < 0.5 ) then BoxSize_new = ( BoxSize ** 3.0 - dV ) ** ( 1.0 / 3.0 ) Increase = 2 dV = -dV else BoxSize_new = ( BoxSize ** 3.0 + dV ) ** ( 1.0 / 3.0 ) Increase = 1 end if dU = 0.0 ZERO = 0.0 CZERO1 = (0.0, 0.0) CZERO2 = (0.0, 0.0) CZERO3 = (0.0, 0.0) CoulCombo = ec * ec * 1.0e10 / ( 4.0 * Pi * eps0 * kB ) BoxRatio = BoxSize_new / BoxSize Xlj_new = Xlj Ylj_new = Ylj Zlj_new = Zlj Xion_new = Xion Yion_new = Yion Zion_new = Zion do j = 1, Nmol lenlj = LENGTHlj( j ) stlj = STARTlj( j ) endlj = stlj + lenlj - 1 lenion = LENGTHion( j ) stion = STARTion( j ) endion = stion + lenion - 1 if( lenion == 0 ) then call Outfold( lenlj, lenion, BoxSize, Xlj_new(stlj:endlj), & Ylj_new(stlj:endlj), Zlj_new(stlj:endlj), & Xion_new, Yion_new, & Zion_new ) else call Outfold( lenlj, lenion, BoxSize, Xlj_new(stlj:endlj), & Ylj_new(stlj:endlj), Zlj_new(stlj:endlj), & Xion_new(stion:endion), Yion_new(stion:endion), & Zion_new(stion:endion) ) end if if( lenion == 0 ) then call CenterOfMass( lenlj, Xlj_new(stlj:endlj), & Ylj_new(stlj:endlj), Zlj_new(stlj:endlj), & TYPElj( stlj:endlj ), Nljgrs, MASSlj, & xcom, ycom, zcom ) else temp1( 1:lenlj ) = Xlj_new( stlj:endlj ) temp1( lenlj+1:lenlj+lenion ) = Xion_new( stion:endion ) temp2( 1:lenlj ) = Ylj_new( stlj:endlj ) temp2( lenlj+1:lenlj+lenion ) = Yion_new( stion:endion ) temp3( 1:lenlj ) = Zlj_new( stlj:endlj ) temp3( lenlj+1:lenlj+lenion ) = Zion_new( stion:endion ) temp4( 1:lenlj ) = TYPElj( stlj:endlj ) temp4( lenlj+1:lenlj+lenion ) = TYPEion( 1:lenion ) + Nljgrs temp5( 1:Nljgrs ) = MASSlj( 1:Nljgrs ) temp5( Nljgrs+1:Nljgrs+Niongrs ) = MASSion( 1:Niongrs ) call CenterOfMass( lenlj + lenion, temp1, temp2, temp3, temp4, & Nljgrs + Niongrs, temp5, xcom, ycom, zcom ) end if Xlj_new( stlj:endlj ) = Xlj_new( stlj:endlj ) + & xcom * ( BoxRatio - 1.0 ) Ylj_new( stlj:endlj ) = Ylj_new( stlj:endlj ) + & ycom * ( BoxRatio - 1.0 ) Zlj_new( stlj:endlj ) = Zlj_new( stlj:endlj ) + & zcom * ( BoxRatio - 1.0 ) do k = stlj, endlj if( Xlj_new(k) > BoxSize_new ) Xlj_new(k) = Xlj_new(k) - & BoxSize_new * aint( Xlj_new(k) / BoxSize_new ) if( Ylj_new(k) > BoxSize_new ) Ylj_new(k) = Ylj_new(k) - & BoxSize_new * aint( Ylj_new(k) / BoxSize_new ) if( Zlj_new(k) > BoxSize_new ) Zlj_new(k) = Zlj_new(k) - & BoxSize_new * aint( Zlj_new(k) / BoxSize_new ) if( Xlj_new(k) < 0.0 ) Xlj_new(k) = Xlj_new(k) - & BoxSize_new * aint( Xlj_new(k) / BoxSize_new - 1 ) if( Ylj_new(k) < 0.0 ) Ylj_new(k) = Ylj_new(k) - & BoxSize_new * aint( Ylj_new(k) / BoxSize_new - 1 ) if( Zlj_new(k) < 0.0 ) Zlj_new(k) = Zlj_new(k) - & BoxSize_new * aint( Zlj_new(k) / BoxSize_new - 1 ) end do if( lenion > 0 ) then Ions = .True. Xion_new( stion:endion ) = Xion_new( stion:endion ) + & xcom * ( BoxRatio - 1.0 ) Yion_new( stion:endion ) = Yion_new( stion:endion ) + & ycom * ( BoxRatio - 1.0 ) Zion_new( stion:endion ) = Zion_new( stion:endion ) + & zcom * ( BoxRatio - 1.0 ) do k = stion, endion if( Xion_new(k) > BoxSize_new ) Xion_new(k) = Xion_new(k) - & BoxSize_new * aint( Xion_new(k) / BoxSize_new ) if( Yion_new(k) > BoxSize_new ) Yion_new(k) = Yion_new(k) - & BoxSize_new * aint( Yion_new(k) / BoxSize_new ) if( Zion_new(k) > BoxSize_new ) Zion_new(k) = Zion_new(k) - & BoxSize_new * aint( Zion_new(k) / BoxSize_new ) if( Xion_new(k) < 0.0 ) Xion_new(k) = Xion_new(k) - & BoxSize_new * aint( Xion_new(k) / BoxSize_new - 1 ) if( Yion_new(k) < 0.0 ) Yion_new(k) = Yion_new(k) - & BoxSize_new * aint( Yion_new(k) / BoxSize_new - 1 ) if( Zion_new(k) < 0.0 ) Zion_new(k) = Zion_new(k) - & BoxSize_new * aint( Zion_new(k) / BoxSize_new - 1 ) end do call SelfMolecule( lenion, Xion_new(stion:endion), Yion_new(stion:endion), & Zion_new(stion:endion), TYPEion(stion:endion), & Nham, Niongrs, CHARGE, BoxSize_new, Alpha, & USELF_MOL_new( j,: ) ) USELF_MOL_new( j,: ) = USELF_MOL_new( j,: ) * CoulCombo else USELF_MOL_new( j,: ) = 0.0 end if end do call e6interact( Nlj, Xlj_new, Ylj_new, Zlj_new, TYPElj, Nmol, LENGTHlj, & Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, BoxSize_new, ULJBOX_new ) NGROUPS = 0 do k = 1, Nlj NGROUPS( TYPElj(k) ) = NGROUPS( TYPElj(k) ) + 1 end do if( CP(1,1,1) > 0.0 ) then EPS_TMP = EPS * CP * ALP / ( ALP - 6.0 ) / 4.0 else EPS_TMP = EPS end if call lrcorr( Nljgrs, Nham, BoxSize_new, NGROUPS, EPS_TMP, SIG, XLRCORR, & ELRCORR, ULJLR_new ) if( Ions ) then call Surf_Move( Nion, Xion_new, Yion_new, Zion_new, TYPEion, & Nham, Niongrs, CHARGE, BoxSize_new, & ZERO, ZERO, ZERO, SUMQX_NEW, SUMQY_NEW, & SUMQZ_NEW, USURF_new ) USURF_new = USURF_new * CoulCombo call Fourier_Move( Nion, Xion_new, Yion_new, Zion_new, & TYPEion, Nham, Niongrs, CHARGE, BoxSize_new, & Kmax, Nkvec, KX, KY, KZ, CONST, & CZERO1, CZERO2, CZERO2, & EXPX_NEW, EXPY_NEW, EXPZ_NEW, & CZERO3, SUMQEXPV_NEW, UFOURIER_new ) UFOURIER_new = UFOURIER_new * CoulCombo call RealInteract( Nion, Xion_new, Yion_new, Zion_new, & TYPEion, Nmol, LENGTHion, Nham, & Niongrs, CHARGE, BoxSize_new, Alpha, UREAL_new ) UREAL_new = UREAL_new * CoulCombo USELF_CH_new = USELF_CH / BoxRatio else USURF_new = 0.0 UFOURIER_new = 0.0 UREAL_new = 0.0 USELF_CH_new = 0.0 end if do h = 1, Nham dU(h,1) = ULJBOX_new(h) + ULJLR_new(h) + USURF_new(h) + & UFOURIER_new(h) + UREAL_new(h) - USELF_CH_new(h) - & sum( USELF_MOL_new( 1:Nmol,h ) ) - & ( ULJBOX(h) + ULJLR(h) + USURF(h) + & UFOURIER(h) + UREAL(h) - USELF_CH(h) - & sum( USELF_MOL( 1:Nmol,h ) ) ) + & Pressure * dV * 1.0e-25 / kB end do TEST_PRESS = exp( real( Nmol) * 3.0 * log( BoxRatio ) - & matmul( BETA, transpose(dU) ) ) dU = - dU BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + Largest PiRatio = LnPiRatio + Nmol * 3.0 * log( BoxRatio ) if( log( ran2(Seed) ) < PiRatio ) then Success = .True. BoxSize = BoxSize_new Xlj = Xlj_new Ylj = Ylj_new Zlj = Zlj_new ULJBOX = ULJBOX_new ULJLR = ULJLR_new LNWSTATE = LNWSTATE + BETA_HAM - LnPiRatio If( Ions ) then Xion = Xion_new Yion = Yion_new Zion = Zion_new USURF = USURF_new UFOURIER = UFOURIER_new UREAL = UREAL_new USELF_CH = USELF_CH_new USELF_MOL = USELF_MOL_new EXPX = EXPX_NEW EXPY = EXPY_NEW EXPZ = EXPZ_NEW SUMQEXPV = SUMQEXPV_NEW SUMQX = SUMQX_NEW SUMQY = SUMQY_NEW SUMQZ = SUMQZ_NEW end if end if return end subroutine Volume_NPT