📄 create.f90
字号:
UINTRA(j) = UINTRA(j) + Uvib + Ubend + Utor
end do
end do
end do
close(30)
else
Nlj = 0
Nion = 0
Nmol_new = 0
SUMQX = 0.0
SUMQY = 0.0
SUMQZ = 0.0
ULJBOX = 0.0
ULJLR = 0.0
ULJ_MOL = 0.0
UINTRA = 0.0
UREAL = 0.0
UFOURIER = 0.0
USELF_CH = 0.0
USURF = 0.0
USELF_MOL = 0.0
UION_MOL = 0.0
EXPX = ( 0.0, 0.0 )
EXPY = ( 0.0, 0.0 )
EXPZ = ( 0.0, 0.0 )
SUMQEXPV = ( 0.0, 0.0 )
PROB_SP(1:Nsp) = real( Nmol(1:Nsp) )
if( Nmol(0) == 0 ) PROB_SP(Nsp) = 1.0
do j = 2, Nsp
PROB_SP(j) = PROB_SP(j-1) + PROB_SP(j)
end do
PROB_SP = PROB_SP / PROB_SP(Nsp)
do while ( Nmol_new(0) < Nmol(0) )
Random = ran2(Seed)
sp = 0
j = 1
do while ( sp == 0 )
if( Random < PROB_SP(j) ) sp = j
j = j + 1
end do
if( Nmol_new(sp) == Nmol(sp) ) cycle
llj = LENLJ(sp)
lion = LENION(sp)
allocate( Xlj_new(llj) )
allocate( Ylj_new(llj) )
allocate( Zlj_new(llj) )
allocate( TYPElj_new(llj) )
allocate( DAMPlj2_new(llj) )
allocate( DAMPlj3_new(llj) )
DAMPlj2_new = 1.0
DAMPlj3_new = 1.0
Count = 0
do j = 1, llj + lion
if( BEADTYPE(j,sp) == 'LJ' ) then
Count = Count + 1
TYPElj_new(Count) = GROUPTYPE(j,sp)
end if
end do
if( lion > 0 ) then
allocate( Xion_new( lion ) )
allocate( Yion_new( lion ) )
allocate( Zion_new( lion ) )
allocate( TYPEion_new( lion ) )
allocate( DAMPion_new( lion ) )
allocate( EXPX_new( 0:Kmax, lion ) )
allocate( EXPY_new( -Kmax:Kmax, lion ) )
allocate( EXPZ_new( -Kmax:Kmax, lion ) )
DAMPion_new = 1.0
Count = 0
do j = 1, llj + lion
if( BEADTYPE(j,sp) == 'ION' ) then
Count = Count + 1
TYPEion_new(Count) = GROUPTYPE(j,sp)
end if
end do
end if
if( Nmol_new(sp) == 0 ) then
BigW_old = 1.0
do j = 1, NSTEPS(sp)
BigW_old = BigW_old * NTRIALS(j,sp)
end do
else
Mol = int( Nmol_new(sp) * ran2(Seed) ) + 1
j = 0
Count = 0
do while ( Count < Mol )
j = j + 1
if( SPECIES(j) == sp ) Count = Count + 1
end do
Mol = j
stlj = STARTlj(Mol)
endlj = stlj + llj - 1
stion = STARTion(Mol)
endion = stion + lion - 1
if( stlj == 1 ) then
Xlj_tmp( 1:Nlj-llj ) = Xlj( endlj+1:Nlj )
Ylj_tmp( 1:Nlj-llj ) = Ylj( endlj+1:Nlj )
Zlj_tmp( 1:Nlj-llj ) = Zlj( endlj+1:Nlj )
TYPElj_tmp( 1:Nlj-llj ) = TYPElj( endlj+1:Nlj )
DAMPlj2_tmp( 1:Nlj-llj ) = DAMPlj2( endlj+1:Nlj )
DAMPlj3_tmp( 1:Nlj-llj ) = DAMPlj3( endlj+1:Nlj )
else
Xlj_tmp( 1:stlj-1 ) = Xlj( 1:stlj-1 )
Xlj_tmp( stlj:Nlj-llj ) = Xlj( endlj+1:Nlj )
Ylj_tmp( 1:stlj-1 ) = Ylj( 1:stlj-1 )
Ylj_tmp( stlj:Nlj-llj ) = Ylj( endlj+1:Nlj )
Zlj_tmp( 1:stlj-1 ) = Zlj( 1:stlj-1 )
Zlj_tmp( stlj:Nlj-llj ) = Zlj( endlj+1:Nlj )
TYPElj_tmp( 1:stlj-1 ) = TYPElj( 1:stlj-1 )
TYPElj_tmp( stlj:Nlj-llj ) = TYPElj( endlj+1:Nlj )
DAMPlj2_tmp( 1:stlj-1 ) = DAMPlj2( 1:stlj-1 )
DAMPlj2_tmp( stlj:Nlj-llj ) = DAMPlj2( endlj+1:Nlj )
DAMPlj3_tmp( 1:stlj-1 ) = DAMPlj3( 1:stlj-1 )
DAMPlj3_tmp( stlj:Nlj-llj ) = DAMPlj3( endlj+1:Nlj )
end if
if( lion > 0 ) then
if( stion == 1 ) then
Xion_tmp( 1:Nion-lion ) = Xion( endion+1:Nion )
Yion_tmp( 1:Nion-lion ) = Yion( endion+1:Nion )
Zion_tmp( 1:Nion-lion ) = Zion( endion+1:Nion )
TYPEion_tmp( 1:Nion-lion ) = TYPEion( endion+1:Nion )
DAMPion_tmp( 1:Nion-lion ) = DAMPion( endion+1:Nion )
else
Xion_tmp( 1:stion-1 ) = Xion( 1:stion-1 )
Xion_tmp( stion:Nion-lion ) = Xion( endion+1:Nion )
Yion_tmp( 1:stion-1 ) = Yion( 1:stion-1 )
Yion_tmp( stion:Nion-lion ) = Yion( endion+1:Nion )
Zion_tmp( 1:stion-1 ) = Zion( 1:stion-1 )
Zion_tmp( stion:Nion-lion ) = Zion( endion+1:Nion )
TYPEion_tmp( 1:stion-1 ) = TYPEion( 1:stion-1 )
TYPEion_tmp( stion:Nion-lion ) = TYPEion( endion+1:Nion )
DAMPion_tmp( 1:stion-1 ) = DAMPion( 1:stion-1 )
DAMPion_tmp( stion:Nion-lion ) = DAMPion( endion+1:Nion )
end if
end if
Xlj_new( 1:llj ) = Xlj( stlj:endlj )
Ylj_new( 1:llj ) = Ylj( stlj:endlj )
Zlj_new( 1:llj ) = Zlj( stlj:endlj )
call e6molecule( llj, Xlj_new, Ylj_new, Zlj_new, &
TYPElj(stlj:endlj), DAMPlj2(stlj:endlj), &
DAMPlj3(stlj:endlj), &
Nlj-llj, Xlj_tmp, Ylj_tmp, Zlj_tmp, &
TYPElj_tmp, DAMPlj2_tmp, DAMPlj3_tmp, &
Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
BoxSize, dULJBOX )
ULJBOX_new = ULJBOX - dULJBOX
NGROUPS = 0
do k = 1, Nlj-llj
NGROUPS( TYPElj_tmp(k) ) = NGROUPS( TYPElj_tmp(k) ) + 1
end do
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, NGROUPS, EPS_CP, SIG, &
XLRCORR, ELRCORR, ULJLR_new )
ULJ_MOL_new = 0.0
Uintra_new = 0.0
if( lion > 0 ) then
Xion_new( 1:lion ) = Xion( stion:endion )
Yion_new( 1:lion ) = Yion( stion:endion )
Zion_new( 1:lion ) = Zion( stion:endion )
call RealMolecule( lion, Xion_new, Yion_new, Zion_new, &
TYPEion(stion:endion), &
DAMPion(stion:endion), &
Nion-lion, Xion_tmp, Yion_tmp, &
Zion_tmp, TYPEion_tmp, DAMPion_tmp, &
Nham, Niongrs, CHARGE, &
BoxSize, Alpha, dUREAL )
UREAL_new = UREAL - dUREAL * CoulCombo
call Surf_Move( lion, -Xion_new, -Yion_new, -Zion_new, &
TYPEion(stion:endion), &
DAMPion(stion:endion), &
Nham, Niongrs, CHARGE, BoxSize, &
SUMQX, SUMQY, SUMQZ, &
SUMQX_new, SUMQY_new, SUMQZ_new, dUSURF )
USURF_new = USURF + dUSURF * CoulCombo
allocate( CZERO1(0:Kmax,lion) )
allocate( CZERO2(-Kmax:Kmax,lion) )
CZERO1 = ( 0.0, 0.0 )
CZERO2 = ( 0.0, 0.0 )
call Fourier_Move( lion, Xion_new, Yion_new, Zion_new, &
TYPEion( stion:endion ), &
DAMPion(stion:endion), &
Nham, Niongrs, -CHARGE, BoxSize, &
Kmax, Nkvec, KX, KY, KZ, CONST, &
CZERO1, CZERO2, CZERO2, &
EXPX_new, EXPY_new, EXPZ_new, &
SUMQEXPV, SUMQEXPV_new, dUFOURIER )
UFOURIER_new = UFOURIER + dUFOURIER * CoulCombo
deallocate( CZERO1 )
deallocate( CZERO2 )
dUSELF_CH = 0.0
do h = 1, Nham
do m = stion, endion
dUSELF_CH(h) = dUSELF_CH(h) + DAMPion(m) * DAMPion(m) * &
CHARGE( TYPEion(m), h ) * &
CHARGE( TYPEion(m), h )
end do
dUSELF_CH(h) = dUSELF_CH(h) * Alpha * CoulCombo / sqrt(Pi) / BoxSize
USELF_CH_new(h) = USELF_CH(h) - dUSELF_CH(h)
end do
USELF_MOL_new = 0.0
UION_MOL_new = 0.0
end if
New = .False.
call Grow( NSTEPS(sp), STEPSTART(:,sp), STEPLENGTH(:,sp), &
NTRIALS(:,sp), llj, lion, MaxBeads, &
METHOD(:,sp), MaxInt, MaxReal, &
INTPARAM(:,:,sp), REALPARAM(:,:,sp), &
BEADTYPE(:,sp), New, 1, &
Xlj_new, Ylj_new, Zlj_new, &
TYPElj(stlj:endlj), DAMPlj2(stlj:endlj), &
DAMPlj3(stlj:endlj), &
Xion_new, Yion_new, Zion_new, &
TYPEion(stion:endion), DAMPion(stion:endion ), &
Nham, ULJBOX_new, ULJLR_new, ULJ_MOL_new, &
UREAL_new, USURF_new, UFOURIER_new, USELF_CH_new, &
USELF_MOL_new, UION_MOL_new, Uintra_new, &
Niongrs, CHARGE, Alpha, &
Kmax, Nkvec, KX, KY, KZ, CONST, &
SUMQX_new, SUMQY_new, SUMQZ_new, &
EXPX_new, EXPY_new, EXPZ_new, &
SUMQEXPV_new, BETA, BIGW_old, LNW, &
Nlj-llj, Xlj_tmp, Ylj_tmp, Zlj_tmp, &
TYPElj_tmp, DAMPlj2_tmp, DAMPlj3_tmp, &
Nion-lion, Xion_tmp, Yion_tmp, Zion_tmp, &
TYPEion_tmp, DAMPion_tmp, &
BoxSize, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
XLRCORR, ELRCORR, &
Nres, Nresmol, reslen, Xresmols, &
Yresmols, Zresmols, Uint_resm, &
Ulj_resm, Uion_resm, &
BONDSAPART(:,:,sp), f14_lj, f14_ion, Seed )
end if
ULJBOX_new = ULJBOX
ULJLR_new = ULJLR
ULJ_MOL_new = 0.0
UREAL_new = UREAL
USURF_new = USURF
UFOURIER_new = UFOURIER
USELF_CH_new = USELF_CH
USELF_MOL_new = 0.0
UION_MOL_new = 0.0
Uintra_new = 0.0
SUMQX_new = SUMQX
SUMQY_new = SUMQY
SUMQZ_new = SUMQZ
SUMQEXPV_new = SUMQEXPV
New = .True.
call Grow( NSTEPS(sp), STEPSTART(:,sp), STEPLENGTH(:,sp), &
NTRIALS(:,sp), llj, lion, MaxBeads, &
METHOD(:,sp), MaxInt, MaxReal, &
INTPARAM(:,:,sp), REALPARAM(:,:,sp), &
BEADTYPE(:,sp), New, 1, &
Xlj_new, Ylj_new, Zlj_new, &
TYPElj_new, DAMPlj2_new, DAMPlj3_new, &
Xion_new, Yion_new, Zion_new, &
TYPEion_new, DAMPion_new, &
Nham, ULJBOX_new, ULJLR_new, ULJ_MOL_new, &
UREAL_new, USURF_new, UFOURIER_new, USELF_CH_new, &
USELF_MOL_new, USELF_MOL_new, Uintra_new, &
Niongrs, CHARGE, Alpha, &
Kmax, Nkvec, KX, KY, KZ, CONST, &
SUMQX_new, SUMQY_new, SUMQZ_new, &
EXPX_new, EXPY_new, EXPZ_new, &
SUMQEXPV_new, BETA, BIGW_new, 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(:,:,sp), f14_lj, f14_ion, Seed )
if( ran2(Seed) < BigW_new(1) / BigW_old(1) ) then ! Approximation
Xlj( Nlj+1:Nlj+llj ) = Xlj_new( 1:llj )
Ylj( Nlj+1:Nlj+llj ) = Ylj_new( 1:llj )
Zlj( Nlj+1:Nlj+llj ) = Zlj_new( 1:llj )
TYPElj( Nlj+1:Nlj+llj ) = TYPElj_new( 1:llj )
DAMPlj2( Nlj+1:Nlj+llj ) = DAMPlj2_new( 1:llj )
DAMPlj3( Nlj+1:Nlj+llj ) = DAMPlj3_new( 1:llj )
ULJBOX = ULJBOX_new
ULJLR = ULJLR_new
ULJ_MOL( Nmol_new(0)+1,: ) = ULJ_MOL_new(:)
if( lion > 0 ) then
Xion( Nion+1:Nion+lion ) = Xion_new( 1:lion )
Yion( Nion+1:Nion+lion ) = Yion_new( 1:lion )
Zion( Nion+1:Nion+lion ) = Zion_new( 1:lion )
TYPEion( Nion+1:Nion+lion ) = TYPEion_new( 1:lion )
DAMPion( Nion+1:Nion+lion ) = DAMPion_new( 1:lion )
UREAL = UREAL_new
USURF = USURF_new
UFOURIER = UFOURIER_new
USELF_CH = USELF_CH_new
SUMQX = SUMQX_new
SUMQY = SUMQY_new
SUMQZ = SUMQZ_new
SUMQEXPV = SUMQEXPV_new
EXPX( :, Nion+1:Nion+lion ) = EXPX_new( :, 1:lion )
EXPY( :, Nion+1:Nion+lion ) = EXPY_new( :, 1:lion )
EXPZ( :, Nion+1:Nion+lion ) = EXPZ_new( :, 1:lion )
USELF_MOL( Nmol_new(0)+1,: ) = USELF_MOL_new(:)
UION_MOL( Nmol_new(0)+1,: ) = UION_MOL_new(:)
end if
UINTRA( Nmol_new(0)+1 ) = Uintra_new
LENGTHlj( Nmol_new(0)+1 ) = llj
LENGTHion( Nmol_new(0)+1 ) = lion
STARTlj( Nmol_new(0)+1 ) = Nlj + 1
STARTion( Nmol_new(0)+1 ) = Nion + 1
SPECIES( Nmol_new(0)+1 ) = sp
Nmol_new(0) = Nmol_new(0) + 1
Nmol_new(sp) = Nmol_new(sp) + 1
Nlj = Nlj + llj
Nion = Nion + lion
end if
deallocate( Xlj_new )
deallocate( Ylj_new )
deallocate( Zlj_new )
deallocate( TYPElj_new )
deallocate( DAMPlj2_new )
deallocate( DAMPlj3_new )
if( lion > 0 ) then
deallocate( Xion_new )
deallocate( Yion_new )
deallocate( Zion_new )
deallocate( TYPEion_new )
deallocate( DAMPion_new )
deallocate( EXPX_new )
deallocate( EXPY_new )
deallocate( EXPZ_new )
end if
end do
end if
do h = 1, Nham
UION(h) = USURF(h) + UFOURIER(h) + UREAL(h) - USELF_CH(h) - &
sum( USELF_MOL( 1:Nmol(0), h) )
ULJ(h) = ULJBOX(h) + ULJLR(h)
UTOT(h) = UION(h) + ULJ(h) + sum( ULJ_MOL( 1:Nmol(0), h) ) + &
sum( UION_MOL( 1:Nmol(0), h) )
end do
LNPSI = 3.0*Nmol(0)*log(BoxSize) - BETA*UTOT
do i = 1, Nsp
LNPSI(:) = LNPSI(:) + Nmol(i)*log(ZETA(i,:))
do k = 1, Nmol(i)
LNPSI = LNPSI - log(real(k))
end do
end do
Largest = maxval( LNW + LNPSI )
LnPi = log( sum( exp( LNW + LNPSI - Largest ) ) ) + Largest
return
end subroutine Create
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -