Subroutine Grow( NSteps, STEPSTART, STEPLENGTH, NTRIALS, lenlj, lenion, & MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, REALPARAM, & BEADTYPE, New, StartGrowthStep, & Xlj_new, Ylj_new, Zlj_new, TYPElj_new, & Xion_new, Yion_new, Zion_new, TYPEion_new, & Nham, ULJBOX, ULJLR, ULJ_MOL, UREAL, USURF, & UFOURIER, USELF_CH, USELF_MOL, UION_MOL, Uintra, & Niongrs, CHARGE, Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, & SUMQX, SUMQY, SUMQZ, EXPX_new, EXPY_new, EXPZ_new, & SUMQEXPV, BigW, NTemp, BETA, LNWSTATE, & Nlj, Xlj, Ylj, Zlj, TYPElj, & Nion, Xion, Yion, Zion, TYPEion, & BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, XLRCORR, ELRCORR, & BONDSAPART, f14_lj, f14_ion, Seed ) implicit none ! This subroutine grows a molecule using configurational bias. ! NSteps is the number of configurational bias steps it takes to grow the molecule. integer, intent(in) :: NSteps ! STEPSTART is the start bead of each CB step. ! STEPLENGTH is the number of beads in each CB step. ! NTRIALS is the number of trial orientations/locations for each CB step. integer, dimension(NSteps) :: STEPSTART integer, dimension(NSteps) :: STEPLENGTH integer, dimension(NSteps) :: NTRIALS ! lenlj is the number of LJ beads in the molecule to be grown. ! lenion is the number of ionic beads in the molecule to be grown. integer, intent(in) :: lenlj, lenion ! METHOD is the method used to grow each bead. ! MaxInt is the maximum number of integer parameters. ! MaxReal is the maximum number of real parameters. ! INTPARAM are the integer parameters needed by the method to grow each bead. ! REALPARAM are the real parameters needed by the method to grow each bead. ! BEADTYPE indicates whether the bead is 'LJ' or 'ION'. integer, intent(in) :: MaxBeads character*10, dimension( MaxBeads ), intent(in) :: METHOD integer, intent(in) :: MaxInt integer, intent(in) :: MaxReal integer, dimension( MaxInt, MaxBeads ), intent(in) :: INTPARAM real, dimension( MaxReal, MaxBeads ), intent(in) :: REALPARAM character*5, dimension( MaxBeads ), intent(in) :: BEADTYPE ! New is a logical to indicate if a new molecule is being grown or the ! Rosenbluth weight of an old molecule is being determined. ! New = .True. when a new molecule is to be grown. logical, intent(in) :: New ! StartGrowthStep is the starting step number for gegrowing the molecule. ! StartGrowthStep = 1, for regrowing the whole molecule. integer, intent(in) :: StartGrowthStep ! Xlj_new, Ylj_new, and Zlj_new contain the new coordinates of the grown ! molecule if New = .True., or the coordinates of an existing molecule if ! New = .False. ! When regrowing part of a molecule Xlj_new, etc. contain the coordinates ! of the molecule up to StartGrowthStep - 1 real, dimension(lenlj), intent(inout) :: Xlj_new real, dimension(lenlj), intent(inout) :: Ylj_new real, dimension(lenlj), intent(inout) :: Zlj_new ! TYPElj_new contains the group identity of the LJ beads to be grown. integer, dimension(lenlj), intent(in) :: TYPElj_new ! Xion_new, Yion_new, and Zion_new have a similar definition to that above for ! Xlj_new, Ylj_new, and Zlj_new. real, dimension(lenion), intent(inout) :: Xion_new real, dimension(lenion), intent(inout) :: Yion_new real, dimension(lenion), intent(inout) :: Zion_new ! TYPEion_new contains the group identity of the ionic beads to be grown. integer, dimension(lenion), intent(in) :: TYPEion_new ! Nham is the number of hamiltonians. integer, intent(in) :: Nham ! ULJBOX is the LJ energy of the phase before and after the molecule is grown. ! ULJLR is the long range correction to the LJ energy. ! ULJ_MOL is the non bonded LJ energy of the molecule to be grown. real, dimension(Nham), intent(inout) :: ULJBOX real, dimension(Nham), intent(inout) :: ULJLR real, dimension(Nham), intent(out) :: ULJ_MOL ! UREAL is the real coulombic energy. ! USURF is the surface energy. ! UFOURIER is the reciprical space coulombic energy. ! USELF_CH is the self charge energy. ! USELF_MOL is the self energy of the molecule to be grown. ! UION_MOL is the non-bonded ionic intramolecular energy of the molecule. real, dimension(Nham), intent(inout) :: UREAL real, dimension(Nham), intent(inout) :: USURF real, dimension(Nham), intent(inout) :: UFOURIER real, dimension(Nham), intent(inout) :: USELF_CH real, dimension(Nham), intent(out) :: USELF_MOL real, dimension(Nham), intent(out) :: UION_MOL ! Uintra is the intramolecular energy of the grown molecule. real, intent(out) :: Uintra ! Niongrs is the number of ion groups. ! CHARGE contains the charge for a given group and hamiltonian. integer, intent(in) :: Niongrs real, dimension(Niongrs, Nham), intent(in) :: CHARGE ! 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 ! SUMQX is the summation of qi * xi for all ions in the box before and after ! the growth. real, dimension(Nham), intent(inout) :: SUMQX, SUMQY, SUMQZ ! EXPX_new, etc. contain exp( i*kx*x ) for the beads of the grown molecule. complex, dimension(0:Kmax, lenion), intent(inout) :: EXPX_new complex, dimension(-Kmax:Kmax, lenion), intent(inout) :: EXPY_new, EXPZ_new ! SUMQEXPV contains the summation of qi*exp(i*(kx*x + ky*y + kz*z)) ! for a given k-vector and hamiltonian before and after the growth. complex, dimension(Nkvec, Nham), intent(inout) :: SUMQEXPV ! BigW is the Rosenbluth factor of the grown molecule. real, intent(out) :: BigW ! NTemp is the number of temperatures being used. integer, intent(in) :: NTemp ! BETA contains the recprical temperatures. real, dimension(NTemp,1), intent(in) :: BETA ! LNWSTATE contains the weight of each state. real, dimension(NTemp, Nham), intent(in) :: LNWSTATE ! Nlj is the number of LJ beads in the phase the molecule is grown in. ! Xlj, Ylj, and Zlj are the coordinates of the Nlj LJ beads. ! TYPElj contains the group identity of the Nlj LJ beads. integer, intent(in) :: Nlj real, dimension(Nlj), intent(in) :: Xlj, Ylj, Zlj integer, dimension(Nlj), intent(in) :: TYPElj ! Nion is the number of ionic beads in the phase the molecule is grown in. ! Xion, Yion, and Zion are the coordinates of the Nion ionic beads. ! TYPEion contains the group identity of the Nion ionic beads. integer, intent(in) :: Nion real, dimension(Nion), intent(in) :: Xion, Yion, Zion integer, dimension(Nion), intent(in) :: TYPEion ! BoxSize is the length of the simulation box. real, intent(in) :: BoxSize ! 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 ! XLRCORR and ELRCORR contain correction factors for the long range LJ energy. real, dimension(605), intent(in) :: XLRCORR real, dimension(605), intent(in) :: ELRCORR ! BONDSAPART gives the bondlengths any two LJ beads are away from each other. integer, dimension(MaxBeads,MaxBeads) :: BONDSAPART ! f14_lj and f14_ion are damping factors for the 1-4 intramolecular interactions real, intent(in) :: f14_lj real, intent(in) :: f14_ion ! Seed is the current random number generator seed value. integer, intent(inout) :: Seed real, external :: ran2 ! Local Variables integer :: h, i, j, k, m integer :: ljb, ionb integer :: ljb_temp, ionb_temp integer :: lenljst, lenionst integer :: Ntry, Nb integer :: Selection integer :: nbendp, noutfold integer :: ntemp1, ntemp2 integer :: ntemp3, ntemp4 integer :: ncount integer, dimension(Nljgrs) :: NLJGROUPS integer, dimension(lenlj+lenion) :: LJBEAD, IONBEAD integer :: nDoOver = 1000 logical :: Success real :: CoulCombo real :: LnPiRatio real :: Largest real :: Random real :: req, rtrial real, dimension(1) :: Zero = 0.0 real, dimension(15) :: Xc, Yc, Zc real, dimension(4) :: ZERO2 real, dimension(Nham) :: Utemp real, dimension(Nham) :: dULJ_MOL_part real, dimension(Nham) :: dUION_MOL_part real, dimension(Nham,1) :: dU real, dimension(NTemp, Nham) :: BETA_HAM real, dimension(NSteps) :: SMALLW real, dimension(Nljgrs,Nljgrs,Nham) :: EPS_TMP real, allocatable, dimension(:) :: PiRatio real, allocatable, dimension(:) :: PROB real, allocatable, dimension(:) :: TempX, TempY, TempZ real, allocatable, dimension(:,:) :: Xlj_tr, Ylj_tr, Zlj_tr real, allocatable, dimension(:,:) :: Xion_tr, Yion_tr, Zion_tr real, allocatable, dimension(:,:) :: SUMQX_tr, SUMQY_tr, SUMQZ_tr real, allocatable, dimension(:,:) :: dULJBOX_tr real, allocatable, dimension(:,:) :: dULJLR_tr real, allocatable, dimension(:,:) :: dULJ_MOL_tr real, allocatable, dimension(:,:) :: dUION_MOL_tr real, allocatable, dimension(:,:) :: dUREAL_tr real, allocatable, dimension(:,:) :: dUSURF_tr real, allocatable, dimension(:,:) :: dUFOURIER_tr real, allocatable, dimension(:,:) :: dUSELF_CH_tr real, allocatable, dimension(:,:) :: dUSELF_MOL_tr real, allocatable, dimension(:,:) :: UVIB real, allocatable, dimension(:,:) :: UBEND real, allocatable, dimension(:,:) :: UTOR real, parameter :: Pi = 3.14159265359 real, parameter :: ec = 1.60217733e-19 real, parameter :: eps0 = 8.854187817e-12 real, parameter :: kB = 1.380658e-23 complex, allocatable, dimension(:,:,:) :: EXPX_tr complex, allocatable, dimension(:,:,:) :: EXPY_tr complex, allocatable, dimension(:,:,:) :: EXPZ_tr complex, allocatable, dimension(:,:,:) :: SUMQEXPV_tr complex, allocatable, dimension(:,:) :: CZERO1, CZERO2 CoulCombo = ec * ec * 1.0e10 / ( 4.0 * Pi * eps0 * kB ) LJBEAD = 0 IONBEAD = 0 ljb = 0 ionb = 0 do i = 1, lenion + lenlj if( BEADTYPE(i) == 'LJ' ) then ljb = ljb + 1 LJBEAD(i) = ljb IONBEAD(i) = ionb else if( BEADTYPE(i) == 'ION' ) then ionb = ionb + 1 IONBEAD(i) = ionb LJBEAD(i) = ljb end if end do ljb = 0 ionb = 0 ! NLJGROUPS, for calculation of ULJLR. NLJGROUPS = 0 ! NLJGROUPS for complete molecules. do i = 1, Nlj NLJGROUPS( TYPElj(i) ) = NLJGROUPS( TYPElj(i) ) + 1 end do ! NLJGROUPS for molecule being grown. do i = 1, StartGrowthStep - 1 do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1 ljb = LJBEAD(k) ionb = IONBEAD(k) if( BEADTYPE(k) == 'LJ' ) then NLJGROUPS( TYPElj_new(ljb) ) = NLJGROUPS( TYPElj_new(ljb) ) + 1 end if end do end do ! Uintra, etc. will continuosly be passed through grow. !Uintra = 0.0 !USELF_MOL = 0.0 !UION_MOL = 0.0 !ULJ_MOL = 0.0 BigW = 1.0 do i = StartGrowthStep, NSteps lenljst = 0 lenionst = 0 do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1 if( BEADTYPE(k) == 'LJ' ) then lenljst = lenljst + 1 else if( BEADTYPE(k) == 'ION' ) then lenionst = lenionst + 1 end if end do Ntry = NTRIALS(i) if( lenljst > 0 ) then allocate( Xlj_tr( Ntry, ljb + 1 : ljb + lenljst ) ) allocate( Ylj_tr( Ntry, ljb + 1 : ljb + lenljst ) ) allocate( Zlj_tr( Ntry, ljb + 1 : ljb + lenljst ) ) end if if( lenionst > 0 ) then allocate( Xion_tr( Ntry, ionb + 1 : ionb + lenionst ) ) allocate( Yion_tr( Ntry, ionb + 1 : ionb + lenionst ) ) allocate( Zion_tr( Ntry, ionb + 1 : ionb + lenionst ) ) allocate( EXPX_tr( Ntry, 0:Kmax, ionb + 1 : ionb + lenionst ) ) allocate( EXPY_tr( Ntry, -Kmax:Kmax, ionb + 1 : ionb + lenionst ) ) allocate( EXPZ_tr( Ntry, -Kmax:Kmax, ionb + 1 : ionb + lenionst ) ) allocate( SUMQEXPV_tr( Ntry, Nkvec, Nham ) ) allocate( SUMQX_tr( Ntry, Nham ) ) allocate( SUMQY_tr( Ntry, Nham ) ) allocate( SUMQZ_tr( Ntry, Nham ) ) end if allocate( dULJBOX_tr( Ntry, Nham ) ) allocate( dULJLR_tr( Ntry, Nham ) ) allocate( dULJ_MOL_tr( Ntry, Nham ) ) allocate( dUREAL_tr( Ntry, Nham ) ) allocate( dUSURF_tr( Ntry, Nham ) ) allocate( dUFOURIER_tr( Ntry, Nham ) ) allocate( dUSELF_CH_tr( Ntry, Nham ) ) allocate( dUSELF_MOL_tr( Ntry, Nham ) ) allocate( dUION_MOL_tr( Ntry, Nham ) ) allocate( UVIB( Ntry, lenlj + lenion ) ) allocate( UBEND( Ntry, lenlj + lenion ) ) allocate( UTOR( Ntry, lenlj + lenion ) ) dULJBOX_tr = 0.0 dULJLR_tr = 0.0 dULJ_MOL_tr = 0.0 dUREAL_tr = 0.0 dUSURF_tr = 0.0 dUFOURIER_tr = 0.0 dUSELF_CH_tr = 0.0 dUSELF_MOL_tr = 0.0 dUION_MOL_tr = 0.0 UVIB = 0.0 UBEND = 0.0 UTOR = 0.0 allocate( PiRatio( Ntry ) ) allocate( PROB( Ntry ) ) SMALLW(i) = 0.0 do j = 1, Ntry ! Find the coordinates of the trial positions. if( .NOT. New .AND. j == 1 ) then do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1 ljb = LJBEAD(k) ionb = IONBEAD(k) if( BEADTYPE(k) == 'LJ' ) then NLJGROUPS( TYPElj_new(ljb) ) = NLJGROUPS( TYPElj_new(ljb) ) + 1 end if if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xlj_new(ljb) Ylj_tr(j,ljb) = Ylj_new(ljb) Zlj_tr(j,ljb) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xion_new(ionb) Yion_tr(j,ionb) = Yion_new(ionb) Zion_tr(j,ionb) = Zion_new(ionb) EXPX_tr(j,:,ionb) = EXPX_new(:,ionb) EXPY_tr(j,:,ionb) = EXPY_new(:,ionb) EXPZ_tr(j,:,ionb) = EXPZ_new(:,ionb) end if select case ( METHOD(k) ) case( 'ConeTor' ) do m = 1, 3 if( BEADTYPE( INTPARAM(m,k) ) == 'LJ' ) then Xc(m) = Xlj_new( LJBEAD( INTPARAM(m,k) ) ) Yc(m) = Ylj_new( LJBEAD( INTPARAM(m,k) ) ) Zc(m) = Zlj_new( LJBEAD( INTPARAM(m,k) ) ) else if( BEADTYPE( INTPARAM(m,k) ) == 'ION' ) then Xc(m) = Xion_new( IONBEAD( INTPARAM(m,k) ) ) Yc(m) = Yion_new( IONBEAD( INTPARAM(m,k) ) ) Zc(m) = Zion_new( IONBEAD( INTPARAM(m,k) ) ) end if end do if( BEADTYPE(k) == 'LJ' ) then Xc(4) = Xlj_new(ljb) Yc(4) = Ylj_new(ljb) Zc(4) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xc(4) = Xion_new(ionb) Yc(4) = Yion_new(ionb) Zc(4) = Zion_new(ionb) end if call Outfold( 4, 0, BoxSize, Xc, Yc, Zc, & ZERO2, ZERO2, ZERO2 ) call IntraTorsion( Xc, Yc, Zc, INTPARAM(4,k), & REALPARAM(3:8,k), UTOR(j,k) ) case( 'Stretch' ) if( BEADTYPE( INTPARAM(1,k) ) == 'LJ' ) then Xc(1) = Xlj_new( LJBEAD( INTPARAM(1,k) ) ) Yc(1) = Ylj_new( LJBEAD( INTPARAM(1,k) ) ) Zc(1) = Zlj_new( LJBEAD( INTPARAM(1,k) ) ) else if( BEADTYPE( INTPARAM(1,k) ) == 'ION' ) then Xc(1) = Xion_new( IONBEAD( INTPARAM(1,k) ) ) Yc(1) = Yion_new( IONBEAD( INTPARAM(1,k) ) ) Zc(1) = Zion_new( IONBEAD( INTPARAM(1,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xc(2) = Xlj_new(ljb) Yc(2) = Ylj_new(ljb) Zc(2) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xc(2) = Xion_new(ionb) Yc(2) = Yion_new(ionb) Zc(2) = Zion_new(ionb) end if call Outfold( 2, 0, BoxSize, Xc, Yc, Zc, & ZERO2, ZERO2, ZERO2 ) rtrial = sqrt( ( Xc(2) - Xc(1) ) ** 2 + ( Yc(2) - Yc(1) ) ** 2 + & ( Zc(2) - Zc(1) ) ** 2 ) UVIB(j,k) = 0.5 * REALPARAM(1,k) * ( rtrial - REALPARAM(2,k) ) ** 2.0 case( 'FxBend' ) if( BEADTYPE( INTPARAM(2,k) ) == 'LJ' ) then Xc(2) = Xlj_new( LJBEAD( INTPARAM(2,k) ) ) Yc(2) = Ylj_new( LJBEAD( INTPARAM(2,k) ) ) Zc(2) = Zlj_new( LJBEAD( INTPARAM(2,k) ) ) else if( BEADTYPE( INTPARAM(2,k) ) == 'ION' ) then Xc(2) = Xion_new( IONBEAD( INTPARAM(2,k) ) ) Yc(2) = Yion_new( IONBEAD( INTPARAM(2,k) ) ) Zc(2) = Zion_new( IONBEAD( INTPARAM(2,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xc(3) = Xlj_new(ljb) Yc(3) = Ylj_new(ljb) Zc(3) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xc(3) = Xion_new(ionb) Yc(3) = Yion_new(ionb) Zc(3) = Zion_new(ionb) end if do m = 1, INTPARAM(1,k) if( BEADTYPE( INTPARAM(m+2,k) ) == 'LJ' ) then Xc(m+3) = Xlj_new( LJBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Ylj_new( LJBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zlj_new( LJBEAD( INTPARAM(m+2,k) ) ) else if( BEADTYPE( INTPARAM(m+2,k) ) == 'ION' ) then Xc(m+3) = Xion_new( IONBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Yion_new( IONBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zion_new( IONBEAD( INTPARAM(m+2,k) ) ) end if end do noutfold = INTPARAM(1,k) + 2 call Outfold( noutfold, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), & ZERO2, ZERO2, ZERO2 ) do m = 1, INTPARAM(1,k) Xc(1) = Xc(m+3) Yc(1) = Yc(m+3) Zc(1) = Zc(m+3) ntemp1 = 2 * m ntemp2 = 1 + 2 * m call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do case( 'FxBendTor' ) if( BEADTYPE( INTPARAM(4,k) ) == 'LJ' ) then Xc(3) = Xlj_new( LJBEAD( INTPARAM(4,k) ) ) Yc(3) = Ylj_new( LJBEAD( INTPARAM(4,k) ) ) Zc(3) = Zlj_new( LJBEAD( INTPARAM(4,k) ) ) else if( BEADTYPE( INTPARAM(4,k) ) == 'ION' ) then Xc(3) = Xion_new( IONBEAD( INTPARAM(4,k) ) ) Yc(3) = Yion_new( IONBEAD( INTPARAM(4,k) ) ) Zc(3) = Zion_new( IONBEAD( INTPARAM(4,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xc(4) = Xlj_new(ljb) Yc(4) = Ylj_new(ljb) Zc(4) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xc(4) = Xion_new(ionb) Yc(4) = Yion_new(ionb) Zc(4) = Zion_new(ionb) end if do m = 1, INTPARAM(1,k) + 2 * INTPARAM(2,k) if( BEADTYPE( INTPARAM(m+4,k) ) == 'LJ' ) then Xc(m+4) = Xlj_new( LJBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Ylj_new( LJBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zlj_new( LJBEAD( INTPARAM(m+4,k) ) ) else if( BEADTYPE( INTPARAM(m+4,k) ) == 'ION' ) then Xc(m+4) = Xion_new( IONBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Yion_new( IONBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zion_new( IONBEAD( INTPARAM(m+4,k) ) ) end if end do noutfold = 2 + INTPARAM(1,k) + 2 * INTPARAM(2,k) call Outfold( noutfold, 0, BoxSize, & Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 ) do m = 1, INTPARAM(1,k) Xc(2) = Xc(m+4) Yc(2) = Yc(m+4) Zc(2) = Zc(m+4) ntemp1 = 2 * m ntemp2 = 1 + 2 * m call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do nbendp = INTPARAM(1,k) do m = 1, INTPARAM(2,k) Xc(1) = Xc(2*m+nbendp+3) Yc(1) = Yc(2*m+nbendp+3) Zc(1) = Zc(2*m+nbendp+3) Xc(2) = Xc(2*m+nbendp+4) Yc(2) = Yc(2*m+nbendp+4) Zc(2) = Zc(2*m+nbendp+4) if( INTPARAM(3,k) == 1 ) then ntemp1 = 2*nbendp-4+6*m ntemp2 = 2*nbendp+1+6*m call IntraTorsion( Xc, Yc, Zc, 1, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) else if( INTPARAM(3,k) == 2 ) then ntemp1 = 2*nbendp-2+4*m ntemp2 = 2*nbendp+1+4*m call IntraTorsion( Xc, Yc, Zc, 2, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) end if UTOR(j,k) = UTOR(j,k) + Utemp(1) end do case( 'StBend' ) if( BEADTYPE( INTPARAM(2,k) ) == 'LJ' ) then Xc(2) = Xlj_new( LJBEAD( INTPARAM(2,k) ) ) Yc(2) = Ylj_new( LJBEAD( INTPARAM(2,k) ) ) Zc(2) = Zlj_new( LJBEAD( INTPARAM(2,k) ) ) else if( BEADTYPE( INTPARAM(2,k) ) == 'ION' ) then Xc(2) = Xion_new( IONBEAD( INTPARAM(2,k) ) ) Yc(2) = Yion_new( IONBEAD( INTPARAM(2,k) ) ) Zc(2) = Zion_new( IONBEAD( INTPARAM(2,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xc(3) = Xlj_new(ljb) Yc(3) = Ylj_new(ljb) Zc(3) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xc(3) = Xion_new(ionb) Yc(3) = Yion_new(ionb) Zc(3) = Zion_new(ionb) end if do m = 1, INTPARAM(1,k) if( BEADTYPE( INTPARAM(m+2,k) ) == 'LJ' ) then Xc(m+3) = Xlj_new( LJBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Ylj_new( LJBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zlj_new( LJBEAD( INTPARAM(m+2,k) ) ) else if( BEADTYPE( INTPARAM(m+2,k) ) == 'ION' ) then Xc(m+3) = Xion_new( IONBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Yion_new( IONBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zion_new( IONBEAD( INTPARAM(m+2,k) ) ) end if end do noutfold = INTPARAM(1,k) + 2 call Outfold( noutfold, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), & ZERO2, ZERO2, ZERO2 ) rtrial = sqrt( ( Xc(3) - Xc(2) ) ** 2 + ( Yc(3) - Yc(2) ) ** 2 + & ( Zc(3) - Zc(2) ) ** 2 ) UVIB(j,k) = 0.5 * REALPARAM(1,k) * ( rtrial - REALPARAM(2,k) ) ** 2.0 do m = 1, INTPARAM(1,k) Xc(1) = Xc(m+3) Yc(1) = Yc(m+3) Zc(1) = Zc(m+3) ntemp1 = 1 + 2 * m ntemp2 = 2 + 2 * m call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do case( 'StBendTor' ) if( BEADTYPE( INTPARAM(4,k) ) == 'LJ' ) then Xc(3) = Xlj_new( LJBEAD( INTPARAM(4,k) ) ) Yc(3) = Ylj_new( LJBEAD( INTPARAM(4,k) ) ) Zc(3) = Zlj_new( LJBEAD( INTPARAM(4,k) ) ) else if( BEADTYPE( INTPARAM(4,k) ) == 'ION' ) then Xc(3) = Xion_new( IONBEAD( INTPARAM(4,k) ) ) Yc(3) = Yion_new( IONBEAD( INTPARAM(4,k) ) ) Zc(3) = Zion_new( IONBEAD( INTPARAM(4,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xc(4) = Xlj_new(ljb) Yc(4) = Ylj_new(ljb) Zc(4) = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then Xc(4) = Xion_new(ionb) Yc(4) = Yion_new(ionb) Zc(4) = Zion_new(ionb) end if do m = 1, INTPARAM(1,k) + 2 * INTPARAM(2,k) if( BEADTYPE( INTPARAM(m+4,k) ) == 'LJ' ) then Xc(m+4) = Xlj_new( LJBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Ylj_new( LJBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zlj_new( LJBEAD( INTPARAM(m+4,k) ) ) else if( BEADTYPE( INTPARAM(m+4,k) ) == 'ION' ) then Xc(m+4) = Xion_new( IONBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Yion_new( IONBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zion_new( IONBEAD( INTPARAM(m+4,k) ) ) end if end do noutfold = 2 + INTPARAM(1,k) + 2 * INTPARAM(2,k) call Outfold( noutfold, 0, BoxSize, & Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 ) rtrial = sqrt( ( Xc(4) - Xc(3) ) ** 2 + ( Yc(4) - Yc(3) ) ** 2 + & ( Zc(4) - Zc(3) ) ** 2 ) UVIB(j,k) = 0.5 * REALPARAM(1,k) * ( rtrial - REALPARAM(2,k) ) ** 2.0 do m = 1, INTPARAM(1,k) Xc(2) = Xc(m+4) Yc(2) = Yc(m+4) Zc(2) = Zc(m+4) ntemp1 = 1 + 2 * m ntemp2 = 2 + 2 * m call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do nbendp = INTPARAM(1,k) do m = 1, INTPARAM(2,k) Xc(1) = Xc(2*m+nbendp+3) Yc(1) = Yc(2*m+nbendp+3) Zc(1) = Zc(2*m+nbendp+3) Xc(2) = Xc(2*m+nbendp+4) Yc(2) = Yc(2*m+nbendp+4) Zc(2) = Zc(2*m+nbendp+4) if( INTPARAM(3,k) == 1 ) then ntemp1 = 2*nbendp-3+6*m ntemp2 = 2*nbendp+2+6*m call IntraTorsion( Xc, Yc, Zc, 1, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) else if( INTPARAM(3,k) == 2 ) then ntemp1 = 2*nbendp-1+4*m ntemp2 = 2*nbendp+2+4*m call IntraTorsion( Xc, Yc, Zc, 2, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) end if UTOR(j,k) = UTOR(j,k) + Utemp(1) end do case default UVIB(j,k) = 0.0 UBEND(j,k) = 0.0 UTOR(j,k) = 0.0 end select end do else if( j == 1 ) then do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1 if( BEADTYPE(k) == 'LJ' ) then NLJGROUPS( TYPElj_new(LJBEAD(k)) ) = NLJGROUPS( TYPElj_new(LJBEAD(k)) ) + 1 end if end do end if k = STEPSTART(i) - 1 cb_step_loop: do while ( k < STEPSTART(i) + STEPLENGTH(i) - 1 ) k = k + 1 ljb = LJBEAD(k) ionb = IONBEAD(k) select case ( METHOD(k) ) case( 'Random' ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = ran2(Seed) * BoxSize Ylj_tr(j,ljb) = ran2(Seed) * BoxSize Zlj_tr(j,ljb) = ran2(Seed) * BoxSize else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = ran2(Seed) * BoxSize Yion_tr(j,ionb) = ran2(Seed) * BoxSize Zion_tr(j,ionb) = ran2(Seed) * BoxSize end if case( 'Sphere' ) if( BEADTYPE( INTPARAM(1,k) ) == 'LJ' ) then Xc(1) = Xlj_new( LJBEAD( INTPARAM(1,k) ) ) Yc(1) = Ylj_new( LJBEAD( INTPARAM(1,k) ) ) Zc(1) = Zlj_new( LJBEAD( INTPARAM(1,k) ) ) else if( BEADTYPE( INTPARAM(1,k) ) == 'ION' ) then Xc(1) = Xion_new( IONBEAD( INTPARAM(1,k) ) ) Yc(1) = Yion_new( IONBEAD( INTPARAM(1,k) ) ) Zc(1) = Zion_new( IONBEAD( INTPARAM(1,k) ) ) end if call Sphere( Xc(1), Yc(1), Zc(1), REALPARAM(1,k), & Xc(2), Yc(2), Zc(2), Seed ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(2) Ylj_tr(j,ljb) = Yc(2) Zlj_tr(j,ljb) = Zc(2) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(2) Yion_tr(j,ionb) = Yc(2) Zion_tr(j,ionb) = Zc(2) end if case( 'Cone' ) do m = 1, 2 if( BEADTYPE( INTPARAM(m,k) ) == 'LJ' ) then Xc(m) = Xlj_new( LJBEAD( INTPARAM(m,k) ) ) Yc(m) = Ylj_new( LJBEAD( INTPARAM(m,k) ) ) Zc(m) = Zlj_new( LJBEAD( INTPARAM(m,k) ) ) else if( BEADTYPE( INTPARAM(m,k) ) == 'ION' ) then Xc(m) = Xion_new( IONBEAD( INTPARAM(m,k) ) ) Yc(m) = Yion_new( IONBEAD( INTPARAM(m,k) ) ) Zc(m) = Zion_new( IONBEAD( INTPARAM(m,k) ) ) end if end do call Outfold( 2, 0, BoxSize, Xc, Yc, Zc, & ZERO2, ZERO2, ZERO2 ) Nb = INTPARAM(3,k) allocate( TempX( Nb ) ) allocate( TempY( Nb ) ) allocate( TempZ( Nb ) ) ntemp1 = Nb+1 ntemp2 = 2*Nb ntemp3 = 2*Nb+1 ntemp4 = 3*Nb call Cone( Nb, Xc, Yc, Zc, REALPARAM(1:Nb,k), & REALPARAM(ntemp1:ntemp2,k), REALPARAM(ntemp3:ntemp4,k), & TempX, TempY, TempZ, Seed ) if( BEADTYPE(k) == 'LJ' ) then ljb_temp = ljb ionb_temp = ionb + 1 else if( BEADTYPE(k) == 'ION' ) then ljb_temp = ljb + 1 ionb_temp = ionb end if do m = 1, Nb if( BEADTYPE( k + m - 1 ) == 'LJ' ) then Xlj_tr(j,ljb_temp) = TempX(m) Ylj_tr(j,ljb_temp) = TempY(m) Zlj_tr(j,ljb_temp) = TempZ(m) ljb_temp = ljb_temp + 1 else if( BEADTYPE( k + m - 1 ) == 'ION' ) then Xion_tr(j,ionb_temp) = TempX(m) Yion_tr(j,ionb_temp) = TempY(m) Zion_tr(j,ionb_temp) = TempZ(m) ionb_temp = ionb_temp + 1 end if end do deallocate( TempX ) deallocate( TempY ) deallocate( TempZ ) case( 'ConeTor' ) do m = 1, 3 if( BEADTYPE( INTPARAM(m,k) ) == 'LJ' ) then Xc(m) = Xlj_new( LJBEAD( INTPARAM(m,k) ) ) Yc(m) = Ylj_new( LJBEAD( INTPARAM(m,k) ) ) Zc(m) = Zlj_new( LJBEAD( INTPARAM(m,k) ) ) else if( BEADTYPE( INTPARAM(m,k) ) == 'ION' ) then Xc(m) = Xion_new( IONBEAD( INTPARAM(m,k) ) ) Yc(m) = Yion_new( IONBEAD( INTPARAM(m,k) ) ) Zc(m) = Zion_new( IONBEAD( INTPARAM(m,k) ) ) end if end do call Outfold( 3, 0, BoxSize, Xc, Yc, Zc, & ZERO2, ZERO2, ZERO2 ) Success = .False. do while( .NOT. Success ) call Cone( 1, Xc(2:3), Yc(2:3), Zc(2:3), & REALPARAM(1,k), Zero, REALPARAM(2,k), & Xc(4), Yc(4), Zc(4), Seed ) call IntraTorsion( Xc, Yc, Zc, INTPARAM(4,k), & REALPARAM(3:8,k), UTOR(j,k) ) dU = - UTOR(j,k) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + & Largest if( log( Ran2(Seed) ) < LnPiRatio ) then Success = .True. if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(4) Ylj_tr(j,ljb) = Yc(4) Zlj_tr(j,ljb) = Zc(4) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(4) Yion_tr(j,ionb) = Yc(4) Zion_tr(j,ionb) = Zc(4) end if end if end do case( 'FxBend' ) if( BEADTYPE( INTPARAM(2,k) ) == 'LJ' ) then Xc(2) = Xlj_new( LJBEAD( INTPARAM(2,k) ) ) Yc(2) = Ylj_new( LJBEAD( INTPARAM(2,k) ) ) Zc(2) = Zlj_new( LJBEAD( INTPARAM(2,k) ) ) else if( BEADTYPE( INTPARAM(2,k) ) == 'ION' ) then Xc(2) = Xion_new( IONBEAD( INTPARAM(2,k) ) ) Yc(2) = Yion_new( IONBEAD( INTPARAM(2,k) ) ) Zc(2) = Zion_new( IONBEAD( INTPARAM(2,k) ) ) end if do m = 1, INTPARAM(1,k) if( BEADTYPE( INTPARAM(m+2,k) ) == 'LJ' ) then Xc(m+3) = Xlj_new( LJBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Ylj_new( LJBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zlj_new( LJBEAD( INTPARAM(m+2,k) ) ) else if( BEADTYPE( INTPARAM(m+2,k) ) == 'ION' ) then Xc(m+3) = Xion_new( IONBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Yion_new( IONBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zion_new( IONBEAD( INTPARAM(m+2,k) ) ) end if end do Xc(3) = Xc(2) Yc(3) = Yc(2) Zc(3) = Zc(2) noutfold = INTPARAM(1,k) + 2 call Outfold( noutfold, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), & ZERO2, ZERO2, ZERO2 ) ! bondlength fixed rtrial = REALPARAM(1,k) Success = .False. ncount = 0 do while( .NOT. Success ) ncount = ncount + 1 if( ncount > nDoOver ) then k = STEPSTART(i) - 1 cycle cb_step_loop end if call Sphere( Xc(2), Yc(2), Zc(2), rtrial, & Xc(3), Yc(3), Zc(3), Seed ) UBEND(j,k) = 0.0 do m = 1, INTPARAM(1,k) Xc(1) = Xc(m+3) Yc(1) = Yc(m+3) Zc(1) = Zc(m+3) ntemp1 = 2 * m ntemp2 = 1 + 2 * m call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do dU = - UBEND(j,k) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + & Largest if( log( Ran2(Seed) ) < LnPiRatio ) then Success = .True. if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(3) Ylj_tr(j,ljb) = Yc(3) Zlj_tr(j,ljb) = Zc(3) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(3) Yion_tr(j,ionb) = Yc(3) Zion_tr(j,ionb) = Zc(3) end if end if end do case( 'FxBendTor' ) if( BEADTYPE( INTPARAM(4,k) ) == 'LJ' ) then Xc(3) = Xlj_new( LJBEAD( INTPARAM(4,k) ) ) Yc(3) = Ylj_new( LJBEAD( INTPARAM(4,k) ) ) Zc(3) = Zlj_new( LJBEAD( INTPARAM(4,k) ) ) else if( BEADTYPE( INTPARAM(4,k) ) == 'ION' ) then Xc(3) = Xion_new( IONBEAD( INTPARAM(4,k) ) ) Yc(3) = Yion_new( IONBEAD( INTPARAM(4,k) ) ) Zc(3) = Zion_new( IONBEAD( INTPARAM(4,k) ) ) end if do m = 1, INTPARAM(1,k) + 2 * INTPARAM(2,k) if( BEADTYPE( INTPARAM(m+4,k) ) == 'LJ' ) then Xc(m+4) = Xlj_new( LJBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Ylj_new( LJBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zlj_new( LJBEAD( INTPARAM(m+4,k) ) ) else if( BEADTYPE( INTPARAM(m+4,k) ) == 'ION' ) then Xc(m+4) = Xion_new( IONBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Yion_new( IONBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zion_new( IONBEAD( INTPARAM(m+4,k) ) ) end if end do Xc(4) = Xc(3) Yc(4) = Yc(3) Zc(4) = Zc(3) noutfold = 2 + INTPARAM(1,k) + 2 * INTPARAM(2,k) call Outfold( noutfold, 0, BoxSize, & Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 ) ! find bondlength min/max rtrial = REALPARAM(1,k) Success = .False. ncount = 0 do while( .NOT. Success ) ncount = ncount + 1 if( ncount > nDoOver ) then k = STEPSTART(i) - 1 cycle cb_step_loop end if call Sphere( Xc(3), Yc(3), Zc(3), rtrial, & Xc(4), Yc(4), Zc(4), Seed ) UBEND(j,k) = 0.0 UTOR(j,k) = 0.0 do m = 1, INTPARAM(1,k) Xc(2) = Xc(m+4) Yc(2) = Yc(m+4) Zc(2) = Zc(m+4) ntemp1 = 2 * m ntemp2 = 1 + 2 * m call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do nbendp = INTPARAM(1,k) do m = 1, INTPARAM(2,k) Xc(1) = Xc(2*m+nbendp+3) Yc(1) = Yc(2*m+nbendp+3) Zc(1) = Zc(2*m+nbendp+3) Xc(2) = Xc(2*m+nbendp+4) Yc(2) = Yc(2*m+nbendp+4) Zc(2) = Zc(2*m+nbendp+4) if( INTPARAM(3,k) == 1 ) then ntemp1 = 2*nbendp-4+6*m ntemp2 = 2*nbendp+1+6*m call IntraTorsion( Xc, Yc, Zc, 1, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) else if( INTPARAM(3,k) == 2 ) then ntemp1 = 2*nbendp-2+4*m ntemp2 = 2*nbendp+1+4*m call IntraTorsion( Xc, Yc, Zc, 2, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) end if UTOR(j,k) = UTOR(j,k) + Utemp(1) end do dU = - UBEND(j,k) - UTOR(j,k) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + & Largest if( log( Ran2(Seed) ) < LnPiRatio ) then success = .True. if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(4) Ylj_tr(j,ljb) = Yc(4) Zlj_tr(j,ljb) = Zc(4) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(4) Yion_tr(j,ionb) = Yc(4) Zion_tr(j,ionb) = Zc(4) end if end if end do case( 'Stretch' ) if( BEADTYPE( INTPARAM(1,k) ) == 'LJ' ) then Xc(1) = Xlj_new( LJBEAD( INTPARAM(1,k) ) ) Yc(1) = Ylj_new( LJBEAD( INTPARAM(1,k) ) ) Zc(1) = Zlj_new( LJBEAD( INTPARAM(1,k) ) ) else if( BEADTYPE( INTPARAM(1,k) ) == 'ION' ) then Xc(1) = Xion_new( IONBEAD( INTPARAM(1,k) ) ) Yc(1) = Yion_new( IONBEAD( INTPARAM(1,k) ) ) Zc(1) = Zion_new( IONBEAD( INTPARAM(1,k) ) ) end if req = REALPARAM(2,k) call bondl( REALPARAM(1,k), req, BETA(1,1), Seed, rtrial ) call Sphere( Xc(1), Yc(1), Zc(1), rtrial, & Xc(2), Yc(2), Zc(2), Seed ) UVIB(j,k) = 0.5 * REALPARAM(1,k) * ( rtrial - req ) ** 2.0 if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(2) Ylj_tr(j,ljb) = Yc(2) Zlj_tr(j,ljb) = Zc(2) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(2) Yion_tr(j,ionb) = Yc(2) Zion_tr(j,ionb) = Zc(2) end if case( 'StBend' ) if( BEADTYPE( INTPARAM(2,k) ) == 'LJ' ) then Xc(2) = Xlj_new( LJBEAD( INTPARAM(2,k) ) ) Yc(2) = Ylj_new( LJBEAD( INTPARAM(2,k) ) ) Zc(2) = Zlj_new( LJBEAD( INTPARAM(2,k) ) ) else if( BEADTYPE( INTPARAM(2,k) ) == 'ION' ) then Xc(2) = Xion_new( IONBEAD( INTPARAM(2,k) ) ) Yc(2) = Yion_new( IONBEAD( INTPARAM(2,k) ) ) Zc(2) = Zion_new( IONBEAD( INTPARAM(2,k) ) ) end if do m = 1, INTPARAM(1,k) if( BEADTYPE( INTPARAM(m+2,k) ) == 'LJ' ) then Xc(m+3) = Xlj_new( LJBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Ylj_new( LJBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zlj_new( LJBEAD( INTPARAM(m+2,k) ) ) else if( BEADTYPE( INTPARAM(m+2,k) ) == 'ION' ) then Xc(m+3) = Xion_new( IONBEAD( INTPARAM(m+2,k) ) ) Yc(m+3) = Yion_new( IONBEAD( INTPARAM(m+2,k) ) ) Zc(m+3) = Zion_new( IONBEAD( INTPARAM(m+2,k) ) ) end if end do Xc(3) = Xc(2) Yc(3) = Yc(2) Zc(3) = Zc(2) noutfold = INTPARAM(1,k) + 2 call Outfold( noutfold, 0, BoxSize, Xc(2:), Yc(2:), Zc(2:), & ZERO2, ZERO2, ZERO2 ) ! find bondlength req = REALPARAM(2,k) call bondl( REALPARAM(1,k), req, BETA(1,1), Seed, rtrial ) UVIB(j,k) = 0.5 * REALPARAM(1,k) * ( rtrial - req ) ** 2.0 Success = .False. ncount = 0 do while( .NOT. Success ) ncount = ncount + 1 if( ncount > nDoOver ) then k = STEPSTART(i) - 1 cycle cb_step_loop end if call Sphere( Xc(2), Yc(2), Zc(2), rtrial, & Xc(3), Yc(3), Zc(3), Seed ) UBEND(j,k) = 0.0 do m = 1, INTPARAM(1,k) Xc(1) = Xc(m+3) Yc(1) = Yc(m+3) Zc(1) = Zc(m+3) ntemp1 = 1 + 2 * m ntemp2 = 2 + 2 * m call IntraBend( Xc, Yc, Zc, REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do dU = - UBEND(j,k) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + & Largest if( log( Ran2(Seed) ) < LnPiRatio ) then Success = .True. if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(3) Ylj_tr(j,ljb) = Yc(3) Zlj_tr(j,ljb) = Zc(3) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(3) Yion_tr(j,ionb) = Yc(3) Zion_tr(j,ionb) = Zc(3) end if end if end do case( 'StBendTor' ) if( BEADTYPE( INTPARAM(4,k) ) == 'LJ' ) then Xc(3) = Xlj_new( LJBEAD( INTPARAM(4,k) ) ) Yc(3) = Ylj_new( LJBEAD( INTPARAM(4,k) ) ) Zc(3) = Zlj_new( LJBEAD( INTPARAM(4,k) ) ) else if( BEADTYPE( INTPARAM(4,k) ) == 'ION' ) then Xc(3) = Xion_new( IONBEAD( INTPARAM(4,k) ) ) Yc(3) = Yion_new( IONBEAD( INTPARAM(4,k) ) ) Zc(3) = Zion_new( IONBEAD( INTPARAM(4,k) ) ) end if do m = 1, INTPARAM(1,k) + 2 * INTPARAM(2,k) if( BEADTYPE( INTPARAM(m+4,k) ) == 'LJ' ) then Xc(m+4) = Xlj_new( LJBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Ylj_new( LJBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zlj_new( LJBEAD( INTPARAM(m+4,k) ) ) else if( BEADTYPE( INTPARAM(m+4,k) ) == 'ION' ) then Xc(m+4) = Xion_new( IONBEAD( INTPARAM(m+4,k) ) ) Yc(m+4) = Yion_new( IONBEAD( INTPARAM(m+4,k) ) ) Zc(m+4) = Zion_new( IONBEAD( INTPARAM(m+4,k) ) ) end if end do Xc(4) = Xc(3) Yc(4) = Yc(3) Zc(4) = Zc(3) noutfold = 2 + INTPARAM(1,k) + 2 * INTPARAM(2,k) call Outfold( noutfold, 0, BoxSize, & Xc(3:), Yc(3:), Zc(3:), ZERO2, ZERO2, ZERO2 ) ! find bondlength req = REALPARAM(2,k) call bondl( REALPARAM(1,k), req, BETA(1,1), Seed, rtrial ) UVIB(j,k) = 0.5 * REALPARAM(1,k) * ( rtrial - req ) ** 2.0 Success = .False. ncount = 0 do while( .NOT. Success ) ncount = ncount + 1 if( ncount > nDoOver ) then k = STEPSTART(i) - 1 cycle cb_step_loop end if call Sphere( Xc(3), Yc(3), Zc(3), rtrial, & Xc(4), Yc(4), Zc(4), Seed ) UBEND(j,k) = 0.0 UTOR(j,k) = 0.0 do m = 1, INTPARAM(1,k) Xc(2) = Xc(m+4) Yc(2) = Yc(m+4) Zc(2) = Zc(m+4) ntemp1 = 1 + 2 * m ntemp2 = 2 + 2 * m call IntraBend( Xc(2:4), Yc(2:4), Zc(2:4), REALPARAM(ntemp1,k), & REALPARAM(ntemp2,k), Utemp(1) ) UBEND(j,k) = UBEND(j,k) + Utemp(1) end do nbendp = INTPARAM(1,k) do m = 1, INTPARAM(2,k) Xc(1) = Xc(2*m+nbendp+3) Yc(1) = Yc(2*m+nbendp+3) Zc(1) = Zc(2*m+nbendp+3) Xc(2) = Xc(2*m+nbendp+4) Yc(2) = Yc(2*m+nbendp+4) Zc(2) = Zc(2*m+nbendp+4) if( INTPARAM(3,k) == 1 ) then ntemp1 = 2*nbendp-3+6*m ntemp2 = 2*nbendp+2+6*m call IntraTorsion( Xc, Yc, Zc, 1, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) else if( INTPARAM(3,k) == 2 ) then ntemp1 = 2*nbendp-1+4*m ntemp2 = 2*nbendp+2+4*m call IntraTorsion( Xc, Yc, Zc, 2, & REALPARAM(ntemp1:ntemp2,k), & Utemp(1) ) end if UTOR(j,k) = UTOR(j,k) + Utemp(1) end do dU = - UBEND(j,k) - UTOR(j,k) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + & Largest if( log( Ran2(Seed) ) < LnPiRatio ) then success = .True. if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(4) Ylj_tr(j,ljb) = Yc(4) Zlj_tr(j,ljb) = Zc(4) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(4) Yion_tr(j,ionb) = Yc(4) Zion_tr(j,ionb) = Zc(4) end if end if end do case( 'Match' ) if( BEADTYPE( INTPARAM(1,k) ) == 'LJ' ) then Xc(1) = Xlj_new( LJBEAD( INTPARAM(1,k) ) ) Yc(1) = Ylj_new( LJBEAD( INTPARAM(1,k) ) ) Zc(1) = Zlj_new( LJBEAD( INTPARAM(1,k) ) ) else if( BEADTYPE( INTPARAM(1,k) ) == 'ION' ) then Xc(1) = Xion_new( IONBEAD( INTPARAM(1,k) ) ) Yc(1) = Yion_new( IONBEAD( INTPARAM(1,k) ) ) Zc(1) = Zion_new( IONBEAD( INTPARAM(1,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xc(1) Ylj_tr(j,ljb) = Yc(1) Zlj_tr(j,ljb) = Zc(1) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xc(1) Yion_tr(j,ionb) = Yc(1) Zion_tr(j,ionb) = Zc(1) end if case( 'Complete' ) end select if( BEADTYPE(k) == 'LJ' ) then if( Xlj_tr(j,ljb) > BoxSize ) Xlj_tr(j,ljb) = Xlj_tr(j,ljb) - & BoxSize * aint( Xlj_tr(j,ljb) / BoxSize ) if( Ylj_tr(j,ljb) > BoxSize ) Ylj_tr(j,ljb) = Ylj_tr(j,ljb) - & BoxSize * aint( Ylj_tr(j,ljb) / BoxSize ) if( Zlj_tr(j,ljb) > BoxSize ) Zlj_tr(j,ljb) = Zlj_tr(j,ljb) - & BoxSize * aint( Zlj_tr(j,ljb) / BoxSize ) if( Xlj_tr(j,ljb) < 0.0 ) Xlj_tr(j,ljb) = Xlj_tr(j,ljb) - & BoxSize * aint( Xlj_tr(j,ljb) / BoxSize - 1 ) if( Ylj_tr(j,ljb) < 0.0 ) Ylj_tr(j,ljb) = Ylj_tr(j,ljb) - & BoxSize * aint( Ylj_tr(j,ljb) / BoxSize - 1 ) if( Zlj_tr(j,ljb) < 0.0 ) Zlj_tr(j,ljb) = Zlj_tr(j,ljb) - & BoxSize * aint( Zlj_tr(j,ljb) / BoxSize - 1 ) Xlj_new(ljb) = Xlj_tr(j,ljb) Ylj_new(ljb) = Ylj_tr(j,ljb) Zlj_new(ljb) = Zlj_tr(j,ljb) else if( BEADTYPE(k) == 'ION' ) then if( Xion_tr(j,ionb) > BoxSize ) Xion_tr(j,ionb) = Xion_tr(j,ionb) - & BoxSize * aint( Xion_tr(j,ionb) / BoxSize ) if( Yion_tr(j,ionb) > BoxSize ) Yion_tr(j,ionb) = Yion_tr(j,ionb) - & BoxSize * aint( Yion_tr(j,ionb) / BoxSize ) if( Zion_tr(j,ionb) > BoxSize ) Zion_tr(j,ionb) = Zion_tr(j,ionb) - & BoxSize * aint( Zion_tr(j,ionb) / BoxSize ) if( Xion_tr(j,ionb) < 0.0 ) Xion_tr(j,ionb) = Xion_tr(j,ionb) - & BoxSize * aint( Xion_tr(j,ionb) / BoxSize - 1 ) if( Yion_tr(j,ionb) < 0.0 ) Yion_tr(j,ionb) = Yion_tr(j,ionb) - & BoxSize * aint( Yion_tr(j,ionb) / BoxSize - 1 ) if( Zion_tr(j,ionb) < 0.0 ) Zion_tr(j,ionb) = Zion_tr(j,ionb) - & BoxSize * aint( Zion_tr(j,ionb) / BoxSize - 1 ) Xion_new(ionb) = Xion_tr(j,ionb) Yion_new(ionb) = Yion_tr(j,ionb) Zion_new(ionb) = Zion_tr(j,ionb) end if end do cb_step_loop end if ! Find the energy of the trial positions. if( lenljst > 0 ) then call e6molecule( lenljst, Xlj_tr(j,:), Ylj_tr(j,:), Zlj_tr(j,:), & TYPElj_new(ljb-lenljst+1:ljb), Nlj, Xlj, Ylj, & Zlj, TYPElj, Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, & BoxSize, dULJBOX_tr(j,:) ) do k = 1, ljb + ionb do m = ljb + ionb - lenljst - lenionst + 1, ljb + ionb if( k < m .AND. BONDSAPART(k,m) >= 3 .AND. & BEADTYPE(k) == 'LJ' .AND. BEADTYPE(m) == 'LJ') then call e6molecule( 1, Xlj_tr(j,LJBEAD(m)), Ylj_tr(j,LJBEAD(m)), & Zlj_tr(j,LJBEAD(m)), TYPElj_new(LJBEAD(m)), & 1, Xlj_new(LJBEAD(k)), Ylj_new(LJBEAD(k)), & Zlj_new(LJBEAD(k)), TYPElj_new(LJBEAD(k)), & Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, & BoxSize, dULJ_MOL_part ) if( BONDSAPART(k,m) == 3 ) dULJ_MOL_part = f14_lj * dULJ_MOL_part dULJ_MOL_tr(j,:) = dULJ_MOL_tr(j,:) + dULJ_MOL_part(:) end if end do 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, NLJGROUPS, EPS_TMP, SIG, & XLRCORR, ELRCORR, Utemp ) dULJLR_tr(j,:) = Utemp(:) - ULJLR(:) end if if( lenionst > 0 ) then call RealMolecule( lenionst, Xion_tr(j,:), Yion_tr(j,:), & Zion_tr(j,:), TYPEion_new(ionb-lenionst+1:ionb), & Nion, Xion, Yion, Zion, TYPEion, Nham, Niongrs, & CHARGE, BoxSize, Alpha, dUREAL_tr(j,:) ) dUREAL_tr(j,:) = dUREAL_tr(j,:) * CoulCombo call Surf_Move( lenionst, Xion_tr(j,:), Yion_tr(j,:), Zion_tr(j,:), & TYPEion_new(ionb-lenionst+1:ionb), Nham, Niongrs, & CHARGE, BoxSize, SUMQX, SUMQY, SUMQZ, SUMQX_tr(j,:), & SUMQY_tr(j,:), SUMQZ_tr(j,:), dUSURF_tr(j,:) ) dUSURF_tr(j,:) = dUSURF_tr(j,:) * CoulCombo allocate( CZERO1(0:Kmax,lenionst) ) allocate( CZERO2(-Kmax:Kmax,lenionst) ) CZERO1 = ( 0.0, 0.0 ) CZERO2 = ( 0.0, 0.0 ) call Fourier_Move( lenionst, Xion_tr(j,:), Yion_tr(j,:), Zion_tr(j,:), & TYPEion_new(ionb-lenionst+1:ionb), Nham, Niongrs, & CHARGE, BoxSize, Kmax, Nkvec, KX, KY, KZ, CONST, & CZERO1, CZERO2, CZERO2, EXPX_tr(j,:,:), & EXPY_tr(j,:,:), EXPZ_tr(j,:,:), SUMQEXPV, & SUMQEXPV_tr(j,:,:), dUFOURIER_tr(j,:) ) dUFOURIER_tr(j,:) = dUFOURIER_tr(j,:) * CoulCombo deallocate( CZERO1 ) deallocate( CZERO2 ) do h = 1, Nham do m = ionb-lenionst+1, ionb dUSELF_CH_tr(j,h) = dUSELF_CH_tr(j,h) + CHARGE( TYPEion_new(m), h ) * & CHARGE( TYPEion_new(m), h ) end do dUSELF_CH_tr(j,h) = dUSELF_CH_tr(j,h) * Alpha * CoulCombo / sqrt(Pi) / BoxSize end do call SelfMolecule( ionb, Xion_new, Yion_new, Zion_new, TYPEion_new, & Nham, Niongrs, CHARGE, BoxSize, Alpha, Utemp ) dUSELF_MOL_tr(j,:) = Utemp(:) * CoulCombo - USELF_MOL(:) do k = 1, ljb + ionb do m = ljb + ionb - lenljst - lenionst + 1, ljb + ionb if( k < m .AND. BONDSAPART(k,m) >= 3 .AND. & BEADTYPE(k) == 'ION' .AND. BEADTYPE(m) == 'ION') then call IonMolecule( 1, Xion_tr(j,IONBEAD(m)), Yion_tr(j,IONBEAD(m)), & Zion_tr(j,IONBEAD(m)), TYPEion_new(IONBEAD(m)), & 1, Xion_new(IONBEAD(k)), Yion_new(IONBEAD(k)), & Zion_new(IONBEAD(k)), TYPEion_new(IONBEAD(k)), & Nham, Niongrs, CHARGE, BoxSize, dUION_MOL_part ) if( BONDSAPART(k,m) == 3 ) dUION_MOL_part = f14_ion * dUION_MOL_part dUION_MOL_tr(j,:) = dUION_MOL_tr(j,:) + dUION_MOL_part(:) * CoulCombo end if end do end do end if dU(:,1) = dULJBOX_tr(j,:) + dULJLR_tr(j,:) + dULJ_MOL_tr(j,:) + & dUREAL_tr(j,:) + dUSURF_tr(j,:) + dUFOURIER_tr(j,:) - & dUSELF_CH_tr(j,:) - dUSELF_MOL_tr(j,:) + dUION_MOL_tr(j,:) dU(:,1) = -dU(:,1) BETA_HAM = matmul( BETA, transpose( dU ) ) Largest = maxval( LNWSTATE + BETA_HAM ) LnPiRatio = log( sum( exp( LNWSTATE + BETA_HAM - Largest ) ) ) + Largest PIRATIO(j) = exp( LnPiRatio ) SMALLW(i) = SMALLW(i) + PIRATIO(j) end do if( SMALLW(i) < 1.0e-200 ) then BigW = 0.0 if( lenljst > 0 ) then deallocate( Xlj_tr ) deallocate( Ylj_tr ) deallocate( Zlj_tr ) end if if( lenionst > 0 ) then deallocate( Xion_tr ) deallocate( Yion_tr ) deallocate( Zion_tr ) deallocate( EXPX_tr ) deallocate( EXPY_tr ) deallocate( EXPZ_tr ) deallocate( SUMQX_tr ) deallocate( SUMQY_tr ) deallocate( SUMQZ_tr ) deallocate( SUMQEXPV_tr ) end if deallocate( dULJBOX_tr ) deallocate( dULJLR_tr ) deallocate( dULJ_MOL_tr ) deallocate( dUREAL_tr ) deallocate( dUSURF_tr ) deallocate( dUFOURIER_tr ) deallocate( dUSELF_CH_tr ) deallocate( dUSELF_MOL_tr ) deallocate( dUION_MOL_tr ) deallocate( UVIB ) deallocate( UBEND ) deallocate( UTOR ) deallocate( PiRatio ) deallocate( PROB ) return end if if( New ) then PROB(1) = PIRATIO(1) / SMALLW(i) do j = 2, Ntry PROB(j) = PROB(j-1) + PIRATIO(j) / SMALLW(i) end do Selection = 0 Random = ran2(Seed) j = 1 do while ( Selection == 0 ) if( Random < PROB(j) ) Selection = j j = j + 1 end do else Selection = 1 end if do k = STEPSTART(i), STEPSTART(i) + STEPLENGTH(i) - 1 ljb = LJBEAD(k) ionb = IONBEAD(k) if( BEADTYPE(k) == 'LJ' ) then Xlj_new(ljb) = Xlj_tr(Selection, ljb) Ylj_new(ljb) = Ylj_tr(Selection, ljb) Zlj_new(ljb) = Zlj_tr(Selection, ljb) else if( BEADTYPE(k) == 'ION' ) then Xion_new(ionb) = Xion_tr(Selection, ionb) Yion_new(ionb) = Yion_tr(Selection, ionb) Zion_new(ionb) = Zion_tr(Selection, ionb) EXPX_new(:,ionb) = EXPX_tr(Selection,:,ionb) EXPY_new(:,ionb) = EXPY_tr(Selection,:,ionb) EXPZ_new(:,ionb) = EXPZ_tr(Selection,:,ionb) end if Uintra = Uintra + UVIB(Selection,k) + UBEND(Selection,k) + UTOR(Selection,k) end do if( lenljst > 0 ) then ULJBOX(:) = ULJBOX(:) + dULJBOX_tr(Selection, :) ULJLR(:) = ULJLR(:) + dULJLR_tr(Selection, :) ULJ_MOL(:) = ULJ_MOL(:) + dULJ_MOL_tr(Selection, :) deallocate( Xlj_tr ) deallocate( Ylj_tr ) deallocate( Zlj_tr ) end if if( lenionst > 0 ) then UREAL(:) = UREAL(:) + dUREAL_tr(Selection, :) USURF(:) = USURF(:) + dUSURF_tr(Selection, :) UFOURIER(:) = UFOURIER(:) + dUFOURIER_tr(Selection, :) USELF_CH(:) = USELF_CH(:) + dUSELF_CH_tr(Selection, :) USELF_MOL(:) = USELF_MOL(:) + dUSELF_MOL_tr(Selection, :) UION_MOL(:) = UION_MOL(:) + dUION_MOL_tr(Selection, :) SUMQX(:) = SUMQX_tr(Selection, :) SUMQY(:) = SUMQY_tr(Selection, :) SUMQZ(:) = SUMQZ_tr(Selection, :) SUMQEXPV(:,:) = SUMQEXPV_tr(Selection,:,:) deallocate( Xion_tr ) deallocate( Yion_tr ) deallocate( Zion_tr ) deallocate( EXPX_tr ) deallocate( EXPY_tr ) deallocate( EXPZ_tr ) deallocate( SUMQX_tr ) deallocate( SUMQY_tr ) deallocate( SUMQZ_tr ) deallocate( SUMQEXPV_tr ) end if BigW = BigW * SMALLW(i) deallocate( dULJBOX_tr ) deallocate( dULJLR_tr ) deallocate( dULJ_MOL_tr ) deallocate( dUREAL_tr ) deallocate( dUSURF_tr ) deallocate( dUFOURIER_tr ) deallocate( dUSELF_CH_tr ) deallocate( dUSELF_MOL_tr ) deallocate( dUION_MOL_tr ) deallocate( UVIB ) deallocate( UBEND ) deallocate( UTOR ) deallocate( PiRatio ) deallocate( PROB ) end do return end Subroutine Grow