📄 grow.f90
字号:
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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -