subroutine Volume_NVT( Nmol, Nlj, Nion, MaxMol, MaxLJ, MaxIon, & Xlj, Ylj, Zlj, TYPElj, Xion, Yion, Zion, TYPEion, & BoxSize, Constant, dvol, & 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, PRESS_PLUS, & PRESS_MINUS, Seed ) implicit none ! This subroutine attempts to change the volume of both simulation boxes ! while keeping the total volume constant. integer, parameter :: NPhases = 2 ! Nmol is the number of molecules in the simulation boxes. ! Nlj is the number of LJ beads in the simulation boxes. ! Nion is the number of ionic beads in the simulation boxes. integer, dimension(NPhases), intent(in) :: Nmol integer, dimension(NPhases), intent(in) :: Nlj integer, dimension(NPhases), intent(in) :: Nion ! MaxMol is the maximum number of molecules per box. ! MaxLJ is the maximum number of LJ beads per box. ! MaxIon is the maximum number of ionic beads per box. integer, intent(in) :: MaxMol integer, intent(in) :: MaxLJ integer, intent(in) :: MaxIon ! Xlj, Ylj, and Zlj are the coordinates of the LJ beads. ! Xion, Yion, and Zion are the coordinates of the ionic beads. real, dimension( MaxLJ, NPhases ), intent(inout) :: Xlj, Ylj, Zlj real, dimension( MaxIon, NPhases ), 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( MaxLJ, NPhases ), intent(in) :: TYPElj integer, dimension( MaxIon, NPhases ), intent(in) :: TYPEion ! BoxSize is the length of the simulation box. real, dimension(NPhases), 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 ! 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(MaxMol, NPhases), intent(in) :: LENGTHlj, LENGTHion integer, dimension(MaxMol, NPhases), 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, MaxIon, NPhases), intent(inout) & :: EXPX complex, dimension(-Kmax:Kmax, MaxIon, NPhases), 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, NPhases), intent(inout) & :: SUMQEXPV ! SUMQX is the summation of qi * xi for all ions in the box. real, dimension(Nham,NPhases), 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,NPhases), intent(inout) :: ULJBOX, ULJLR real, dimension(Nham,NPhases), intent(inout) :: UREAL, UFOURIER real, dimension(Nham,NPhases), intent(inout) :: USURF, USELF_CH real, dimension(MaxMol,Nham,NPhases), intent(inout) :: USELF_MOL ! Increase is a flag to indicate whether the volume of box 1 or 2 was Increased. ! Increase = 1 for Box 1. ! Increase = 2 for Box 2. integer, intent(out) :: Increase ! Success is a logical indicating whether the move was successful or not. logical, intent(out) :: Success real, dimension(NTemp, Nham), intent(out) :: PRESS_PLUS, PRESS_MINUS ! 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(MaxLJ + MaxIon) :: temp4 integer :: h, i, j, k logical :: Ions = .False. real :: CoulCombo real :: Largest, dV real :: PiRatio, LnPiRatio real :: xcom, ycom, zcom real, dimension(NPhases) :: BoxSize_new, BoxRatio real, dimension(MaxLJ, NPhases) :: Xlj_new, Ylj_new, Zlj_new real, dimension(MaxLJ + MaxIon) :: temp1, temp2, temp3 real, dimension(Nljgrs + Niongrs) :: temp5 real, dimension(MaxIon, NPhases) :: Xion_new, Yion_new, Zion_new real, dimension(NTemp, Nham) :: BETA_HAM real, dimension(Nham,NPhases) :: ULJBOX_new, ULJLR_new real, dimension(Nham,NPhases) :: UREAL_new, UFOURIER_new real, dimension(Nham,NPhases) :: USURF_new, USELF_CH_new real, dimension(Nham,NPhases) :: dUPHASE real, dimension(MaxMol,Nham,NPhases) :: USELF_MOL_new real, dimension(Nham,NPhases) :: SUMQX_NEW, SUMQY_NEW, SUMQZ_NEW real, dimension(Nham,NPhases) :: 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 , MaxIon , NPhases) :: EXPX_NEW complex, dimension(-Kmax:Kmax , MaxIon , NPhases) :: EXPY_NEW complex, dimension(-Kmax:Kmax , MaxIon , NPhases) :: EXPZ_NEW complex, dimension(Nkvec, Nham, NPhases) :: SUMQEXPV_NEW complex, dimension(Nkvec, Nham) :: CZERO3 complex, allocatable, dimension(:,:) :: CZERO1 complex, allocatable, dimension(:,:) :: CZERO2 Success = .False. CZERO3 = ( 0.0, 0.0 ) dV = dvol if( .NOT. Constant ) then dV = dV * ran2(Seed) * ( BoxSize(1) ** 3.0 + BoxSize(2) ** 3.0 ) end if if( ran2(Seed) < 0.5 ) then BoxSize_new(1) = ( BoxSize(1) ** 3.0 - dV ) ** ( 1.0 / 3.0 ) BoxSize_new(2) = ( BoxSize(2) ** 3.0 + dV ) ** ( 1.0 / 3.0 ) Increase = 2 else BoxSize_new(1) = ( BoxSize(1) ** 3.0 + dV ) ** ( 1.0 / 3.0 ) BoxSize_new(2) = ( BoxSize(2) ** 3.0 - dV ) ** ( 1.0 / 3.0 ) Increase = 1 end if dU = 0.0 dUPHASE = 0.0 PRESS_PLUS = 1.0 PRESS_MINUS = 1.0 ZERO = 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 i = 1, NPhases do j = 1, Nmol(i) lenlj = LENGTHlj( j, i ) stlj = STARTlj( j, i ) endlj = stlj + lenlj - 1 lenion = LENGTHion( j, i ) stion = STARTion( j, i ) endion = stion + lenion - 1 call Outfold( lenlj, lenion, BoxSize(i), Xlj_new(stlj:endlj,i), & Ylj_new(stlj:endlj,i), Zlj_new(stlj:endlj,i), & Xion_new(stion:endion,i), Yion_new(stion:endion,i), & Zion_new(stion:endion,i) ) if( lenion == 0 ) then call CenterOfMass( lenlj, Xlj_new(stlj:endlj,i), & Ylj_new(stlj:endlj,i), Zlj_new(stlj:endlj,i), & TYPElj( stlj:endlj,i ), Nljgrs, MASSlj, & xcom, ycom, zcom ) else temp1( 1:lenlj ) = Xlj_new( stlj:endlj, i ) temp1( lenlj+1:lenlj+lenion ) = Xion_new( stion:endion, i ) temp2( 1:lenlj ) = Ylj_new( stlj:endlj, i ) temp2( lenlj+1:lenlj+lenion ) = Yion_new( stion:endion, i ) temp3( 1:lenlj ) = Zlj_new( stlj:endlj, i ) temp3( lenlj+1:lenlj+lenion ) = Zion_new( stion:endion, i ) temp4( 1:lenlj ) = TYPElj( stlj:endlj, i ) temp4( lenlj+1:lenlj+lenion ) = TYPEion( 1:lenion, i ) + 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, i ) = Xlj_new( stlj:endlj, i ) + & xcom * ( BoxRatio(i) - 1.0 ) Ylj_new( stlj:endlj, i ) = Ylj_new( stlj:endlj, i ) + & ycom * ( BoxRatio(i) - 1.0 ) Zlj_new( stlj:endlj, i ) = Zlj_new( stlj:endlj, i ) + & zcom * ( BoxRatio(i) - 1.0 ) do k = stlj, endlj if( Xlj_new(k,i) > BoxSize_new(i) ) Xlj_new(k,i) = Xlj_new(k,i) - & BoxSize_new(i) * aint( Xlj_new(k,i) / BoxSize_new(i) ) if( Ylj_new(k,i) > BoxSize_new(i) ) Ylj_new(k,i) = Ylj_new(k,i) - & BoxSize_new(i) * aint( Ylj_new(k,i) / BoxSize_new(i) ) if( Zlj_new(k,i) > BoxSize_new(i) ) Zlj_new(k,i) = Zlj_new(k,i) - & BoxSize_new(i) * aint( Zlj_new(k,i) / BoxSize_new(i) ) if( Xlj_new(k,i) < 0.0 ) Xlj_new(k,i) = Xlj_new(k,i) - & BoxSize_new(i) * aint( Xlj_new(k,i) / BoxSize_new(i) - 1 ) if( Ylj_new(k,i) < 0.0 ) Ylj_new(k,i) = Ylj_new(k,i) - & BoxSize_new(i) * aint( Ylj_new(k,i) / BoxSize_new(i) - 1 ) if( Zlj_new(k,i) < 0.0 ) Zlj_new(k,i) = Zlj_new(k,i) - & BoxSize_new(i) * aint( Zlj_new(k,i) / BoxSize_new(i) - 1 ) end do if( lenion > 0 ) then Ions = .True. Xion_new( stion:endion, i ) = Xion_new( stion:endion, i ) + & xcom * ( BoxRatio(i) - 1.0 ) Yion_new( stion:endion, i ) = Yion_new( stion:endion, i ) + & ycom * ( BoxRatio(i) - 1.0 ) Zion_new( stion:endion, i ) = Zion_new( stion:endion, i ) + & zcom * ( BoxRatio(i) - 1.0 ) do k = stion, endion if( Xion_new(k,i) > BoxSize_new(i) ) Xion_new(k,i) = Xion_new(k,i) - & BoxSize_new(i) * aint( Xion_new(k,i) / BoxSize_new(i) ) if( Yion_new(k,i) > BoxSize_new(i) ) Yion_new(k,i) = Yion_new(k,i) - & BoxSize_new(i) * aint( Yion_new(k,i) / BoxSize_new(i) ) if( Zion_new(k,i) > BoxSize_new(i) ) Zion_new(k,i) = Zion_new(k,i) - & BoxSize_new(i) * aint( Zion_new(k,i) / BoxSize_new(i) ) if( Xion_new(k,i) < 0.0 ) Xion_new(k,i) = Xion_new(k,i) - & BoxSize_new(i) * aint( Xion_new(k,i) / BoxSize_new(i) - 1 ) if( Yion_new(k,i) < 0.0 ) Yion_new(k,i) = Yion_new(k,i) - & BoxSize_new(i) * aint( Yion_new(k,i) / BoxSize_new(i) - 1 ) if( Zion_new(k,i) < 0.0 ) Zion_new(k,i) = Zion_new(k,i) - & BoxSize_new(i) * aint( Zion_new(k,i) / BoxSize_new(i) - 1 ) end do call SelfMolecule( lenion, Xion_new(stion:endion,i), Yion_new(stion:endion,i), & Zion_new(stion:endion,i), TYPEion(stion:endion,i), & Nham, Niongrs, CHARGE, BoxSize_new(i), Alpha, & USELF_MOL_new( j,:,i ) ) USELF_MOL_new( j,:,i ) = USELF_MOL_new( j,:,i ) * CoulCombo else USELF_MOL_new( j,:,i ) = 0.0 end if end do call e6interact( Nlj(i), Xlj_new(:,i), Ylj_new(:,i), Zlj_new(:,i), & TYPElj(:,i), Nmol(i), LENGTHlj(:,i), & Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, & BoxSize_new(i), ULJBOX_new(:,i) ) NGROUPS = 0 do k = 1, Nlj(i) NGROUPS( TYPElj(k,i) ) = NGROUPS( TYPElj(k,i) ) + 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(i), NGROUPS, EPS_TMP, SIG, XLRCORR, & ELRCORR, ULJLR_new(:,i) ) if( Ions ) then call Surf_Move( Nion(i), Xion_new(:,i), Yion_new(:,i), Zion_new(:,i), & TYPEion(:,i), Nham, Niongrs, CHARGE, BoxSize_new(i), & ZERO, ZERO, ZERO, SUMQX_NEW(:, i), SUMQY_NEW(:, i), & SUMQZ_NEW(:, i), USURF_new(:,i) ) USURF_new(:,i) = USURF_new(:,i) * CoulCombo allocate( CZERO1( 0:Kmax, Nion(i) ) ) allocate( CZERO2( -Kmax:Kmax, Nion(i) ) ) CZERO1 = ( 0.0, 0.0 ) CZERO2 = ( 0.0, 0.0 ) call Fourier_Move( Nion(i), Xion_new(:,i), Yion_new(:,i), Zion_new(:,i), & TYPEion(:,i), Nham, Niongrs, CHARGE, BoxSize_new(i), & Kmax, Nkvec, KX, KY, KZ, CONST, CZERO1, CZERO2, CZERO2, & EXPX_NEW(:,1:Nion(i),i), EXPY_NEW(:,1:Nion(i),i), & EXPZ_NEW(:,1:Nion(i),i), CZERO3, & SUMQEXPV_NEW(:,:,i), UFOURIER_new(:,i) ) UFOURIER_new(:,i) = UFOURIER_new(:,i) * CoulCombo deallocate( CZERO1 ) deallocate( CZERO2 ) call RealInteract( Nion(i), Xion_new(:,i), Yion_new(:,i), Zion_new(:,i), & TYPEion(:,i), Nmol(i), LENGTHion(:,i), Nham, & Niongrs, CHARGE, BoxSize_new(i), Alpha, UREAL_new(:,i) ) UREAL_new(:,i) = UREAL_new(:,i) * CoulCombo USELF_CH_new(:,i) = USELF_CH(:,i) / BoxRatio(i) else USURF_new(:,i) = 0.0 UFOURIER_new(:,i) = 0.0 UREAL_new(:,i) = 0.0 USELF_CH_new(:,i) = 0.0 end if do h = 1, Nham dU(h,1) = ULJBOX_new(h,i) + ULJLR_new(h,i) + USURF_new(h,i) + & UFOURIER_new(h,i) + UREAL_new(h,i) - USELF_CH_new(h,i) - & sum( USELF_MOL_new( 1:Nmol(i),h,i ) ) - & ( ULJBOX(h,i) + ULJLR(h,i) + USURF(h,i) + & UFOURIER(h,i) + UREAL(h,i) - USELF_CH(h,i) - & sum( USELF_MOL( 1:Nmol(i),h,i ) ) ) dUPHASE(h,i) = dU(h,1) end do if( Increase == i ) then PRESS_PLUS = exp( real( Nmol(i) ) * 3.0 * log( BoxRatio(i) ) - & matmul( BETA, transpose(dU) ) ) else PRESS_MINUS = exp( real( Nmol(i) ) * 3.0 * log( BoxRatio(i) ) - & matmul( BETA, transpose(dU) ) ) end if end do dU(:,1) = - ( dUPHASE(:, 1) + dUPHASE(:, 2) ) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + Largest PiRatio = LnPiRatio + Nmol(1) * 3.0 * log( BoxRatio(1) ) + & Nmol(2) * 3.0 * log( BoxRatio(2) ) 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_NVT