📄 alpha_ch.f90
字号:
else if( stlj + lenlj - 1 == Nlj ) then
call e6molecule2( lenlj, Xlj_old, Ylj_old, Zlj_old, TYPElj(stlj:endlj), &
DAMPlj2_old, DAMPlj3_old, DAMPlj2_new, DAMPlj3_new, &
Nlj - lenlj, Xlj(1:stlj-1), Ylj(1:stlj-1), &
Zlj(1:stlj-1), TYPElj(1:stlj-1), &
DAMPlj2(1:stlj-1), DAMPlj3(1:stlj-1), &
DAMPlj2(1:stlj-1), DAMPlj3(1:stlj-1), &
Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
BoxSize, dULJBOX )
else
temp1( 1:stlj-1 ) = Xlj( 1:stlj-1 )
temp1( stlj:Nlj-lenlj ) = Xlj( endlj+1:Nlj )
temp2( 1:stlj-1 ) = Ylj( 1:stlj-1)
temp2( stlj:Nlj-lenlj ) = Ylj( endlj+1:Nlj)
temp3( 1:stlj-1 ) = Zlj( 1:stlj-1)
temp3( stlj:Nlj-lenlj ) = Zlj( endlj+1:Nlj)
temp4( 1:stlj-1 ) = TYPElj( 1:stlj-1)
temp4( stlj:Nlj-lenlj ) = TYPElj( endlj+1:Nlj)
temp9( 1:stlj-1 ) = DAMPlj2( 1:stlj-1)
temp9( stlj:Nlj-lenlj ) = DAMPlj2( endlj+1:Nlj)
temp10( 1:stlj-1 ) = DAMPlj3( 1:stlj-1)
temp10( stlj:Nlj-lenlj ) = DAMPlj3( endlj+1:Nlj)
call e6molecule2( lenlj, Xlj_old, Ylj_old, Zlj_old, TYPElj(stlj:endlj), &
DAMPlj2_old, DAMPlj3_old, DAMPlj2_new, DAMPlj3_new, &
Nlj - lenlj, temp1, temp2, temp3, temp4, &
temp9, temp10, temp9, temp10, &
Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
BoxSize, dULJBOX )
end if
! Calculation of ULJ_MOL.
dULJ_MOL(:) = 0.0
! 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 grow.f90.
do k = 1, lenlj + lenion - 1
do m = k + 1, lenlj + lenion
if( BONDSAPART(k,m,SpID) >= 3 .AND. BEADTYPE(k,SpID) == 'LJ' .AND. &
BEADTYPE(m,SpID) == 'LJ') then
ljbk = LJBEAD(k)
ljbm = LJBEAD(m)
call e6molecule( 1, Xlj_old(ljbm), Ylj_old(ljbm), &
Zlj_old(ljbm), TYPElj(stlj+ljbm-1), &
DAMPlj2_new(ljbm), DAMPlj3_new(ljbm), &
1, Xlj_old(ljbk), Ylj_old(ljbk), &
Zlj_old(ljbk), TYPElj(stlj+ljbk-1), &
DAMPlj2_new(ljbk), DAMPlj3_new(ljbk), &
Nham, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
BoxSize, U_part )
if( BONDSAPART(k,m,SpID) == 3 ) U_part = f14_lj * U_part
dULJ_MOL = dULJ_MOL + U_part
end if
end do
end do
dULJ_MOL(:) = dULJ_MOL(:) - ULJ_MOL(tag_mol,:)
! Stop commenting lines here.
! Coulombic Energies
if( lenion > 0 ) then
dDAMPion = DAMPion_new - DAMPion_old
max_damp = maxval( abs(dDAMPion) )
end if
if( lenion > 0 .AND. max_damp > 1.0e-7 ) then
CoulCombo = ec * ec * 1.0e10 / ( 4.0 * Pi * eps0 * kB )
if( stion == 1 ) then
if( Nmol(0) == 1 ) then
dUREAL = 0.0
else
call RealMolecule2( lenion, Xion_old, Yion_old, Zion_old, &
TYPEion(1:lenion), &
DAMPion_old, DAMPion_new, &
Nion - lenion, Xion(endion+1:Nion), &
Yion(endion+1:Nion), Zion(endion+1:Nion), &
TYPEion(endion+1:Nion), &
DAMPion(endion+1:Nion), DAMPion(endion+1:Nion), &
Nham, Niongrs, CHARGE, BoxSize, Alpha, dUREAL )
end if
else if( stion + lenion - 1 == Nion ) then
call RealMolecule2( lenion, Xion_old, Yion_old, Zion_old, &
TYPEion(stion:endion), &
DAMPion_old, DAMPion_new, &
Nion - lenion, Xion(1:stion-1), &
Yion(1:stion-1), Zion(1:stion-1), &
TYPEion(1:stion-1), &
DAMPion(1:stion-1), DAMPion(1:stion-1), &
Nham, Niongrs, CHARGE, BoxSize, Alpha, dUREAL )
else
temp5( 1:stion-1 ) = Xion( 1:stion-1 )
temp5( stion:Nion-lenion ) = Xion( endion+1:Nion )
temp6( 1:stion-1 ) = Yion( 1:stion-1 )
temp6( stion:Nion-lenion ) = Yion( endion+1:Nion )
temp7( 1:stion-1 ) = Zion( 1:stion-1 )
temp7( stion:Nion-lenion ) = Zion( endion+1:Nion )
temp8( 1:stion-1 ) = TYPEion( 1:stion-1 )
temp8( stion:Nion-lenion ) = TYPEion( endion+1:Nion )
temp11( 1:stion-1 ) = DAMPion( 1:stion-1 )
temp11( stion:Nion-lenion ) = DAMPion( endion+1:Nion )
call RealMolecule2( lenion, Xion_old, Yion_old, Zion_old, &
TYPEion(stion:endion), &
DAMPion_old, DAMPion_new, &
Nion - lenion, temp5, temp6, temp7, &
temp8, temp11, temp11, &
Nham, Niongrs, CHARGE, BoxSize, Alpha, dUREAL )
end if
dUREAL = dUREAL * CoulCombo
call Fourier_Move2( lenion, TYPEion( stion:endion ), &
dDAMPion, &
Nham, Niongrs, CHARGE, BoxSize, &
Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX(:,stion:endion), &
EXPY(:,stion:endion), &
EXPZ(:,stion:endion), &
SUMQEXPV, SUMQEXPV_NEW, &
dUFOURIER )
dUFOURIER = dUFOURIER * CoulCombo
call Surf_Move( lenion, Xion_old, Yion_old, Zion_old, &
TYPEion( stion:endion ), dDAMPion, &
Nham, Niongrs, CHARGE, BoxSize, &
SUMQX, SUMQY, SUMQZ, SUMQX_NEW, &
SUMQY_NEW, SUMQZ_NEW, dUSURF )
dUSURF = dUSURF * CoulCombo
call SelfMolecule( lenion, Xion_old, Yion_old, Zion_old, &
TYPEion(stion:endion), DAMPion_new, &
Nham, Niongrs, CHARGE, BoxSize, Alpha, &
dUSELF_MOL )
dUSELF_MOL(:) = dUSELF_MOL(:) * CoulCombo
dUSELF_MOL = dUSELF_MOL - USELF_MOL(tag_mol,:)
dUION_MOL = 0.0
! 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 grow.f90.
do k = 1, lenlj + lenion - 1
do m = k + 1, lenlj + lenion
if( BONDSAPART(k,m,SpID) >= 3 .AND. BEADTYPE(k,SpID) == 'ION' .AND. &
BEADTYPE(m,SpID) == 'ION') then
ionbk = IONBEAD(k)
ionbm = IONBEAD(m)
call IonMolecule( 1, Xion_old(ionbm), Yion_old(ionbm), &
Zion_old(ionbm), TYPEion(stion+ionbm-1), &
DAMPion_new(ionbm), &
1, Xion_old(ionbk), Yion_old(ionbk), &
Zion_old(ionbk), TYPEion(stion+ionbk-1), &
DAMPion_new(ionbk), &
Nham, Niongrs, CHARGE, BoxSize, U_part )
if( BONDSAPART(k,m,SpID) == 3 ) U_part = f14_ion * U_part
dUION_MOL = dUION_MOL + U_part
end if
end do
end do
dUION_MOL = dUION_MOL * CoulCombo
dUION_MOL(:) = dUION_MOL(:) - UION_MOL(tag_mol,:)
! Stop commenting lines here.
dUSELF_CH = 0.0
do h = 1, Nham
do j = stion, endion
dUSELF_CH(h) = dUSELF_CH(h) + &
CHARGE( TYPEion(j), h ) * CHARGE( TYPEion(j), h ) * &
( DAMPion_new(j-stion+1) * DAMPion_new(j-stion+1) - &
DAMPion_old(j-stion+1) * DAMPion_old(j-stion+1) )
end do
dUSELF_CH(h) = dUSELF_CH(h) * &
Alpha / sqrt(Pi) / BoxSize * CoulCombo
end do
else
dUREAL = 0.0
dUSURF = 0.0
dUFOURIER = 0.0
dUSELF_MOL = 0.0
dUION_MOL = 0.0
dUSELF_CH = 0.0
end if
dU = dULJBOX + dULJ_MOL + &
dUREAL + dUSURF + dUFOURIER - dUSELF_CH - dUSELF_MOL + dUION_MOL
ZETA_tmp(:) = ZETA(SpID,:)
LNPSI_new = LNPSI - BETA * dU
Largest = maxval( LNW + LNPSI_new )
LnPi_new = log( sum( exp( LNW + LNPSI_new - Largest ) ) ) + Largest
if( updown == 1 ) then
LnPi_tmp = LnPi_new + ( EEW(tag_val + 1, SpID) - EEW(tag_val, SpID) )
else
LnPi_tmp = LnPi_new - ( EEW(tag_val, SpID) - EEW(tag_val - 1, SpID) )
end if
if( log( ran2(Seed) ) < LnPi_tmp - LnPi ) then
Success = .true.
LnPi = LnPi_new
LNPSI = LNPSI_new
ULJBOX = ULJBOX + dULJBOX
ULJ_MOL(tag_mol,:) = ULJ_MOL(tag_mol,:) + dULJ_MOL(:)
DAMPlj2( stlj:endlj ) = DAMPlj2_new( 1:lenlj )
DAMPlj3( stlj:endlj ) = DAMPlj3_new( 1:lenlj )
if( lenion > 0 ) then
SUMQEXPV = SUMQEXPV_NEW
SUMQX = SUMQX_NEW
SUMQY = SUMQY_NEW
SUMQZ = SUMQZ_NEW
UFOURIER = UFOURIER + dUFOURIER
UREAL = UREAL + dUREAL
USURF = USURF + dUSURF
USELF_CH = USELF_CH + dUSELF_CH
USELF_MOL(tag_mol,:) = USELF_MOL(tag_mol,:) + dUSELF_MOL(:)
UION_MOL(tag_mol,:) = UION_MOL(tag_mol,:) + dUION_MOL(:)
DAMPion( stion:endion ) = DAMPion_new( 1:lenion )
end if
tag_val = tag_val + updown
if( tag_val == EESTEPS(SpID) ) then
tag_val = 0
tag_mol = 0
end if
else
if( tag_val == EESTEPS(SpID) ) then
tag_val = 0
tag_mol = 0
end if
end if
deallocate( Xlj_old )
deallocate( Ylj_old )
deallocate( Zlj_old )
deallocate( DAMPlj2_old )
deallocate( DAMPlj3_old )
deallocate( DAMPlj2_new )
deallocate( DAMPlj3_new )
if( lenion > 0 ) then
deallocate( Xion_old )
deallocate( Yion_old )
deallocate( Zion_old )
deallocate( DAMPion_old )
deallocate( DAMPion_new )
deallocate( dDAMPion )
end if
return
end subroutine alpha_ch
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -