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, DAMPlj2_new, DAMPlj3_new, & Xion_new, Yion_new, Zion_new, & TYPEion_new, DAMPion_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, BETA, BIGW, LNW, & Nlj, Xlj, Ylj, Zlj, & TYPElj, DAMPlj2, DAMPlj3, & Nion, Xion, Yion, Zion, & TYPEion, DAMPion, & BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, & XLRCORR, ELRCORR, & Nres, Nresmol, reslen, Xresmols, & Yresmols, Zresmols, Uint_resm, & Ulj_resm, Uion_resm, & 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 real, dimension(lenlj), intent(inout) :: DAMPlj2_new real, dimension(lenlj), intent(inout) :: DAMPlj3_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 real, dimension(lenion), intent(inout) :: DAMPion_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(inout) :: 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. 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(inout) :: USELF_MOL real, dimension(Nham), intent(inout) :: UION_MOL ! Uintra is the intramolecular energy of the grown molecule. real, intent(inout) :: 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 ! beta contains the recprical temperatures. real, dimension(Nham), intent(in) :: BETA ! BIGW is the Rosenbluth factor of the grown molecule. real, dimension(Nham), intent(out) :: BIGW ! LNW contains the weight of each state. real, dimension(Nham), intent(in) :: LNW ! 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 real, dimension(Nlj), intent(in) :: DAMPlj2, DAMPlj3 ! 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 real, dimension(Nion), intent(in) :: DAMPion ! 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 is a rank 3 array containing the C_ij parameters for each hamiltonian. ! ALP is a rank 3 array containing the alpha_ij parameters for each hamiltonian. ! RMAX is a rank 3 array containing 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 integer, intent(in) :: Nres, Nresmol integer, dimension(Nres), intent(in) :: reslen real, dimension(Nresmol, MaxBeads, Nres), intent(in) :: Xresmols, Yresmols, Zresmols real, dimension(Nresmol, Nres), intent(in) :: Uint_resm real, dimension(Nresmol, Nham, Nres), intent(in) :: Ulj_resm, Uion_resm ! 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, n integer :: ljb, ionb integer :: ljb_temp, ionb_temp integer :: lenljst, lenionst integer :: Ntry, Nb integer :: Selection integer :: ncount integer :: FxSt, newold integer :: nDoOver = 1000 integer :: resl, resstart integer :: resmol integer, dimension(Nljgrs) :: NLJGROUPS integer, dimension(lenlj+lenion) :: LJBEAD, IONBEAD real :: CoulCombo real :: Random real :: xtr, ytr, ztr real, dimension(Nham) :: Utemp real, dimension(Nham) :: dULJ_MOL_part real, dimension(Nham) :: dUION_MOL_part real, dimension(Nham) :: dU real, dimension(Nham) :: W real, dimension(Nham,NSteps) :: SMALLW real, dimension(Nljgrs, Nljgrs, Nham) :: EPS_CP real, dimension(MaxBeads) :: Xres, Yres, Zres real, allocatable, dimension(:,:) :: BOLTZ 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(:,:) :: Xres_tr, Yres_tr, Zres_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(:,:) :: 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(:,:) :: dUION_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 ) W = exp( LNW ) 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 resstart = 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( Xres_tr( Ntry, MaxBeads ) ) allocate( Yres_tr( Ntry, MaxBeads ) ) allocate( Zres_tr( Ntry, MaxBeads ) ) 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( BOLTZ( Nham, 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 if( resstart > 0 ) then do m = 1, resl Xres_tr(j,m) = Xres(m) Yres_tr(j,m) = Yres(m) Zres_tr(j,m) = Zres(m) end do end if select case ( METHOD(k) ) case( 'Stretch' ) FxSt = 2 newold = 2 if( BEADTYPE(k) == 'LJ' ) then xtr = Xlj_new(ljb) ytr = Ylj_new(ljb) ztr = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then xtr = Xion_new(ionb) ytr = Yion_new(ionb) ztr = Zion_new(ionb) end if call Stretch( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, & FxSt, newold, & xtr, ytr, ztr, & UVIB(j,k), Seed ) case( 'FxBend' ) FxSt = 1 newold = 2 if( BEADTYPE(k) == 'LJ' ) then xtr = Xlj_new(ljb) ytr = Ylj_new(ljb) ztr = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then xtr = Xion_new(ionb) ytr = Yion_new(ionb) ztr = Zion_new(ionb) end if call Bend( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), Seed ) case( 'FxBendTor' ) FxSt = 1 newold = 2 if( BEADTYPE(k) == 'LJ' ) then xtr = Xlj_new(ljb) ytr = Ylj_new(ljb) ztr = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then xtr = Xion_new(ionb) ytr = Yion_new(ionb) ztr = Zion_new(ionb) end if call BendTor( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), UTOR(j,k), & Seed ) case( 'StBend' ) FxSt = 2 newold = 2 if( BEADTYPE(k) == 'LJ' ) then xtr = Xlj_new(ljb) ytr = Ylj_new(ljb) ztr = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then xtr = Xion_new(ionb) ytr = Yion_new(ionb) ztr = Zion_new(ionb) end if call Bend( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), Seed ) case( 'StBendTor' ) FxSt = 2 newold = 2 if( BEADTYPE(k) == 'LJ' ) then xtr = Xlj_new(ljb) ytr = Ylj_new(ljb) ztr = Zlj_new(ljb) else if( BEADTYPE(k) == 'ION' ) then xtr = Xion_new(ionb) ytr = Yion_new(ionb) ztr = Zion_new(ionb) end if call BendTor( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), UTOR(j,k), & Seed ) case( 'ResRand' ) n = INTPARAM(1,k) if( resstart == 0 ) resstart = k resl = reslen(n) do m = 1, resl if( BEADTYPE(m+resstart-1) == 'LJ' ) then Xres(m) = Xlj_new( LJBEAD(m+resstart-1) ) Yres(m) = Ylj_new( LJBEAD(m+resstart-1) ) Zres(m) = Zlj_new( LJBEAD(m+resstart-1) ) else if( BEADTYPE(m+resstart-1) == 'ION' ) then Xres(m) = Xion_new( IONBEAD(m+resstart-1) ) Yres(m) = Yion_new( IONBEAD(m+resstart-1) ) Zres(m) = Zion_new( IONBEAD(m+resstart-1) ) end if end do do m = 1, resl Xres_tr(j,m) = Xres(m) Yres_tr(j,m) = Yres(m) Zres_tr(j,m) = Zres(m) 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' ) FxSt = 1 newold = 1 call Stretch( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, & FxSt, newold, & xtr, ytr, ztr, & UVIB(j,k), Seed ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if case( 'Cone' ) Nb = INTPARAM(3,k) allocate( TempX( Nb ) ) allocate( TempY( Nb ) ) allocate( TempZ( Nb ) ) call cone2( lenlj, lenion, Nb, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, & 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( 'FxBend' ) FxSt = 1 newold = 1 call Bend( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), Seed ) if( ncount /= 0 ) then k = STEPSTART(i) - 1 cycle cb_step_loop else if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if end if case( 'FxBendTor' ) FxSt = 1 newold = 1 call BendTor( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), UTOR(j,k), & Seed ) if( ncount /= 0 ) then k = STEPSTART(i) - 1 cycle cb_step_loop else if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if end if case( 'Stretch' ) FxSt = 2 newold = 1 call Stretch( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, & FxSt, newold, & xtr, ytr, ztr, & UVIB(j,k), Seed ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if case( 'StBend' ) FxSt = 2 newold = 1 call Bend( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), Seed ) if( ncount /= 0 ) then k = STEPSTART(i) - 1 cycle cb_step_loop else if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if end if case( 'StBendTor' ) FxSt = 2 newold = 1 call BendTor( lenlj, lenion, & Xlj_new, Ylj_new, Zlj_new, & Xion_new, Yion_new, Zion_new, & LJBEAD, IONBEAD, & BEADTYPE, MaxInt, MaxReal, & INTPARAM(:,k), REALPARAM(:,k), & BoxSize, Nham, BETA, LNW, & FxSt, newold, & ncount, nDoOver, & xtr, ytr, ztr, & UVIB(j,k), UBEND(j,k), UTOR(j,k), & Seed ) if( ncount /= 0 ) then k = STEPSTART(i) - 1 cycle cb_step_loop else if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if end if case( 'Match' ) if( BEADTYPE( INTPARAM(1,k) ) == 'LJ' ) then xtr = Xlj_new( LJBEAD( INTPARAM(1,k) ) ) ytr = Ylj_new( LJBEAD( INTPARAM(1,k) ) ) ztr = Zlj_new( LJBEAD( INTPARAM(1,k) ) ) else if( BEADTYPE( INTPARAM(1,k) ) == 'ION' ) then xtr = Xion_new( IONBEAD( INTPARAM(1,k) ) ) ytr = Yion_new( IONBEAD( INTPARAM(1,k) ) ) ztr = Zion_new( IONBEAD( INTPARAM(1,k) ) ) end if if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = xtr Ylj_tr(j,ljb) = ytr Zlj_tr(j,ljb) = ztr else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = xtr Yion_tr(j,ionb) = ytr Zion_tr(j,ionb) = ztr end if case( 'Complete' ) case( 'ResRand' ) n = INTPARAM(1,k) if( j == 1 ) then if( resstart == 0 ) resstart = k resl = reslen(n) resmol = int( ran2(Seed) * Nresmol ) + 1 do m = 1, resl Xres(m) = Xresmols(resmol,m,n) Yres(m) = Yresmols(resmol,m,n) Zres(m) = Zresmols(resmol,m,n) end do end if if(new) UTOR(j,k) = Uint_resm(resmol,n) if(new) dULJ_MOL_tr(j,:) = Ulj_resm(resmol,:,n) if(new) dUION_MOL_tr(j,:) = Uion_resm(resmol,:,n) call ResRand( resl, Xres, Yres, Zres, & BoxSize, Seed ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xres(1) Ylj_tr(j,ljb) = Yres(1) Zlj_tr(j,ljb) = Zres(1) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xres(1) Yion_tr(j,ionb) = Yres(1) Zion_tr(j,ionb) = Zres(1) end if do m = 1, resl Xres_tr(j,m) = Xres(m) Yres_tr(j,m) = Yres(m) Zres_tr(j,m) = Zres(m) end do case( 'ResRot' ) call ResRot( resl, Xres, Yres, Zres, & BoxSize, Seed ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xres(2) Ylj_tr(j,ljb) = Yres(2) Zlj_tr(j,ljb) = Zres(2) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xres(2) Yion_tr(j,ionb) = Yres(2) Zion_tr(j,ionb) = Zres(2) end if do m = 1, resl Xres_tr(j,m) = Xres(m) Yres_tr(j,m) = Yres(m) Zres_tr(j,m) = Zres(m) end do case( 'ResCir' ) call ResCir( resl, Xres, Yres, Zres, & BoxSize, Seed ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xres(3) Ylj_tr(j,ljb) = Yres(3) Zlj_tr(j,ljb) = Zres(3) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xres(3) Yion_tr(j,ionb) = Yres(3) Zion_tr(j,ionb) = Zres(3) end if do m = 1, resl Xres_tr(j,m) = Xres(m) Yres_tr(j,m) = Yres(m) Zres_tr(j,m) = Zres(m) end do case( 'ResMatch' ) if( BEADTYPE(k) == 'LJ' ) then Xlj_tr(j,ljb) = Xres_tr(j,k-resstart+1) Ylj_tr(j,ljb) = Yres_tr(j,k-resstart+1) Zlj_tr(j,ljb) = Zres_tr(j,k-resstart+1) else if( BEADTYPE(k) == 'ION' ) then Xion_tr(j,ionb) = Xres_tr(j,k-resstart+1) Yion_tr(j,ionb) = Yres_tr(j,k-resstart+1) Zion_tr(j,ionb) = Zres_tr(j,k-resstart+1) end if 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 if( resstart > 0 ) then do m = 1, resl if( Xres_tr(j,m) > BoxSize ) Xres_tr(j,m) = Xres_tr(j,m) - & BoxSize * aint( Xres_tr(j,m) / BoxSize ) if( Yres_tr(j,m) > BoxSize ) Yres_tr(j,m) = Yres_tr(j,m) - & BoxSize * aint( Yres_tr(j,m) / BoxSize ) if( Zres_tr(j,m) > BoxSize ) Zres_tr(j,m) = Zres_tr(j,m) - & BoxSize * aint( Zres_tr(j,m) / BoxSize ) if( Xres_tr(j,m) < 0.0 ) Xres_tr(j,m) = Xres_tr(j,m) - & BoxSize * aint( Xres_tr(j,m) / BoxSize - 1 ) if( Yres_tr(j,m) < 0.0 ) Yres_tr(j,m) = Yres_tr(j,m) - & BoxSize * aint( Yres_tr(j,m) / BoxSize - 1 ) if( Zres_tr(j,m) < 0.0 ) Zres_tr(j,m) = Zres_tr(j,m) - & BoxSize * aint( Zres_tr(j,m) / BoxSize - 1 ) Xres(m) = Xres_tr(j,m) Yres(m) = Yres_tr(j,m) Zres(m) = Zres_tr(j,m) end do 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), & DAMPlj2_new(ljb-lenljst+1:ljb), & DAMPlj3_new(ljb-lenljst+1:ljb), & Nlj, Xlj, Ylj, Zlj, & TYPElj, DAMPlj2, DAMPlj3, & Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, & BoxSize, dULJBOX_tr(j,:) ) ! If the intramolecular nonbonded interactions are included ! as part of the reservoir energy, then comment the ! following lines. Also change the appropriate lines in ! resread.f90, cranksh.f90 and alpha_ch.f90. 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)), & DAMPlj2_new(LJBEAD(m)), DAMPlj3_new(LJBEAD(m)), & 1, Xlj_new(LJBEAD(k)), Ylj_new(LJBEAD(k)), & Zlj_new(LJBEAD(k)), TYPElj_new(LJBEAD(k)), & DAMPlj2_new(LJBEAD(k)), DAMPlj3_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 ! Stop commenting lines here. do h = 1, Nham if( CP(1,1,h) > 0.0 ) then EPS_CP(:,:,h) = EPS(:,:,h) * CP(:,:,h) * ALP(:,:,h) / & ( ALP(:,:,h) - 6.0 ) / 4.0 else EPS_CP(:,:,h) = EPS(:,:,h) * ALP(:,:,h) / ( ALP(:,:,h) - 6.0 ) * & ( ALP(:,:,h) / 6.0 ) ** ( 6.0 / ( ALP(:,:,h) - 6.0 ) ) / 4.0 end if end do call lrcorr( Nljgrs, Nham, BoxSize, NLJGROUPS, EPS_CP, 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), & DAMPion_new(ionb-lenionst+1:ionb), & Nion, Xion, Yion, Zion, TYPEion, DAMPion, & 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), & DAMPion_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), & DAMPion_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) + & DAMPion_new(m) * & DAMPion_new(m) * & 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, DAMPion_new, & Nham, Niongrs, CHARGE, BoxSize, Alpha, Utemp ) dUSELF_MOL_tr(j,:) = Utemp(:) * CoulCombo - USELF_MOL(:) ! If the intramolecular nonbonded interactions are included ! as part of the reservoir energy, then comment the ! following lines. Also change the appropriate lines in ! resread.f90, cranksh.f90 and alpha_ch.f90. 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)), & DAMPion_new(IONBEAD(m)), & 1, Xion_new(IONBEAD(k)), Yion_new(IONBEAD(k)), & Zion_new(IONBEAD(k)), TYPEion_new(IONBEAD(k)), & DAMPion_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 ! Stop commenting lines here. end if dU(:) = 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 = -dU * BETA BOLTZ(:,j) = exp( dU(:) ) SMALLW(:,i) = SMALLW(:,i) + BOLTZ(:,j) end do ! End of loop over number of trial locations, j if( minval( SMALLW( 1:Nham,i ) ) < 1.0e-250 ) 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( Xres_tr ) deallocate( Yres_tr ) deallocate( Zres_tr ) 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( BOLTZ ) deallocate( PROB ) return end if if( New ) then PROB(1) = sum( BOLTZ( 1:Nham,1 ) * W( 1:Nham ) ) / & sum( SMALLW( 1:Nham,i ) * W( 1:Nham ) ) do j = 2, Ntry PROB(j) = PROB(j-1) + sum( BOLTZ( 1:Nham,j ) * W( 1:Nham ) ) / & sum( SMALLW( 1:Nham,i ) * W( 1:Nham ) ) 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 do m = 1, resl Xres(m) = Xres_tr(selection,m) Yres(m) = Yres_tr(selection,m) Zres(m) = Zres_tr(selection,m) end do deallocate( Xres_tr ) deallocate( Yres_tr ) deallocate( Zres_tr ) 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( BOLTZ ) deallocate( PROB ) end do return end Subroutine Grow