📄 gcmc.f90
字号:
RES_ATT, RES_SUC, &
BEADTYPE(:,i), GROUPTYPE(:,i), &
BONDSAPART(:,:,i), Nljgrs, &
EPS, SIG, CP, ALP, RMAX, &
Niongrs, CHARGE, &
f14_lj, f14_ion, &
BETA, LNW, BoxSize, &
UResInt(i), UResLj(:,i), UResIon(:,i), &
Seed, th_h, phi_h )
end do
if( Nres > 0 ) then
write(*,*) RES_ATT(1), RES_SUC(1)
write(*,*) RES_ATT(2), RES_SUC(2)
write(*,*) ' Uintra = ', sum(Uint_resm) / real(Nresmol)
write(*,*) ' Ulj = ', sum(Ulj_resm) / real(Nresmol)
write(*,*) ' Uion = ', sum(Uion_resm) / real(Nresmol)
write(*,*) ' Time = ', timef() - tstart, ' sec.' ! # C_ALPHA
!write(*,*) ' Time = ', etime(xxx) - tstart, ' sec.' ! # U_ALPHA
end if
open( 50, file = PlotsFile )
close(50, status = 'delete' )
open( 90, file = WFile )
close(90, status = 'delete' )
if( histtype == 'list' ) then
UN_HistFile = trim(UN_HistFile)//'.dat'
do i = 1, MaxSp
N_HistFile(i) = trim(N_HistFile(0))//'_'//char(48+i)//'.dat'
end do
allocate( Npart( Nstorage/Ncollect + 1, MaxSp ) )
allocate( Energy( Nstorage/Ncollect + 1 ) )
else
N_HistFile(1) = trim(N_HistFile(0))//'.dat'
end if
if( histtype == 'matrix' ) then
inquire (file = N_HistFile(1), exist = exists)
if( .NOT. exists ) then
Nsim_min = MaxMol
Nsim_max = 0
end if
if( FromDisk .AND. Nequil > 0 ) then
Urange = Umax - Umin
Urange = Urange / ( 1.0 + 2.0 * SafetyFactor )
Usim_min = Umin + SafetyFactor * Urange
Usim_max = Umax - SafetyFactor * Urange
Ubins = 0
else if( FromDisk .AND. Nequil == 0 .AND. exists ) then
Usim_min = 1.0e10
Usim_max = -1.0e10
allocate( UN_HIST( Ubins, 0:MaxMol, Nham ) )
UN_HIST = 0.0
open(80, file = N_HistFile(1))
do i = Nsim_min(1), Nsim_max(1)
read(80,*) j, N_HIST(i,1,1:Nham)
end do
close(80)
do h = 1, Nham
filename = trim(UN_HistFile)//'_h'//char( 48 + h/10 )//char(48 + mod(h,10))//'.dat'
open( 70, file = filename )
read(70,*)
read(70,*)
do i = Nsim_min(1), Nsim_max(1)
read(70,*) j, j, estart
k = nint( ( estart - 0.5*Uwidth - Umin ) / Uwidth ) + 1
read(70,*) UN_HIST( k:k+j-1, i, h )
end do
close(70)
end do
else if( FromDisk .AND. Nequil == 0 .AND. .NOT. exists ) then
allocate( UN_HIST( Ubins, 0:MaxMol, Nham ) )
UN_HIST = 0.0
end if
end if
New = .True.
call banner( ResultsFile, ' ', Infostring, Hostname, Datestring, &
Timestring, New )
call WriteData( InitialConf, Nham, Nljgrs, Niongrs, &
Nsp, MaxSp, MaxSteps, MaxBeads, MaxInt, MaxReal, &
BETA, ZETA, NAMElj, MASSlj, EPS, SIG, CP, ALP, KAPPA, &
LAMDA, NAMEion, MASSion, CHARGE, NAMEsp, NSTEPS, &
LENLJ, LENION, BEADTYPE, GROUPTYPE, NTRIALS, &
STEPSTART, STEPLENGTH, METHOD, INTPARAM, REALPARAM, &
ERSTEPS, ERSTART, EREND, EESTEPS, MaxEE, BEADDAMP, &
FromDisk, Nmol, Volume, nmoves, PROB_MOVE, PROB_DR, &
PROB_SP_CD, PROB_SP_RG, NinCycle, Nstorage, Nadjust, &
Ncollect, Uwidth, Seed, ResultsFile )
open( 10, file = 'lrcorr.dat' )
do i = 1, 605
read(10,*) XLRCORR(i), ELRCORR(i)
end do
close(10)
if( histtype == 'list' ) then
open(100, file=UN_HistFile )
write(100, "(1x,3(G11.5), 3(F10.4,2x))") &
1.0/BETA(1), log(ZETA(1:Nsp,1)) / BETA(1), &
BoxSize, BoxSize, BoxSize
close(100)
Nsim_max = 0
Nsim_min = MaxMol
end if
call Fourier_Setup( Kmax, Alpha, CONST, KX, KY, KZ, Nkvec )
allocate( SUMQEXPV( Nkvec, Nham ) )
call Create( InitialConf, Nham, Nljgrs, Niongrs, Nsp, &
MaxSp, MaxSteps, MaxBeads, &
MaxInt, MaxReal, MaxEE, &
BETA, ZETA, EPS, SIG, CP, ALP, RMAX, CHARGE, &
NSTEPS, LENLJ, LENION, BEADTYPE, GROUPTYPE, &
BEADDAMP, NTRIALS, STEPSTART, STEPLENGTH, &
METHOD, INTPARAM , REALPARAM , BONDSAPART, &
EESTEPS, FromDisk, BoxSize, &
Nlj, Nion, MaxMol, MaxLJ, MaxIon, Nmol, &
Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, TYPEion, DAMPion, &
SPECIES, LENGTHlj, LENGTHion, &
STARTlj, STARTion, &
LNW, LNPSI, LnPi, XLRCORR, ELRCORR, &
Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, &
SUMQX, SUMQY, SUMQZ, &
ULJBOX, ULJLR, UREAL, UFOURIER, &
USURF, USELF_CH, USELF_MOL, &
ULJ_MOL, UION_MOL, UINTRA, &
tag_val, tag_mol, &
Nres, Nresmol, reslen, Xresmols, &
Yresmols, Zresmols, Uint_resm, &
Ulj_resm, Uion_resm, &
f14_lj, f14_ion, Seed )
write(*,*)
write(*,"(A)") ' The system has been created. '
NATT = 0
NSUC = 0
TOT_NATT = 0
TOT_NSUC = 0
do h = 1, Nham
UION(h) = USURF(h) + UFOURIER(h) + UREAL(h) - USELF_CH(h) - &
sum( USELF_MOL( 1:Nmol(0), h) ) + &
sum( UION_MOL( 1:Nmol(0), h) )
ULJ(h) = ULJBOX(h) + ULJLR(h) + sum( ULJ_MOL( 1:Nmol(0), h) )
UTOT(h) = UION(h) + ULJ(h)
end do
write(*,*)
write(*,"(A,T15,I8)") ' MC Step = ', StartConf
write(*,"(A,T15,F15.4)") ' Rho = ', &
dot_product( Nmol(1:Nsp) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / Volume
write(*,11) ' UTotal = ', UTOT(1:Nham)
write(*,11) ' LNPSI = ', LNPSI(1:Nham)
write(*,"(A,T15,4I8)") ' Nmol = ', Nmol(1:Nsp)
write(*,"(A,T15,4I8)") ' Tag Value = ', tag_val(1:Nsp)
write(*,"(A,T15,4I8)") ' Tag Mol. = ', tag_mol(1:Nsp)
11 format( A, T15, <Nham>F15.4 )
Istore = 1
nt1 = Nmol(0)
do i = StartConf + 1, StartConf + Nequil + Nprod
Moves = Moves + 1
Ranq = ran2(Seed)
if( Ranq < PROB_MOVE(1) ) then
tt = timef() ! # C_ALPHA
! tt = etime(xxx) ! # U_ALPHA
if( Nmol(0) > 0 ) then
PROB_SP_DR(1) = Nmol(1)
do k = 2, Nsp
PROB_SP_DR(k) = PROB_SP_DR(k-1) + Nmol(k)
end do
PROB_SP_DR = PROB_SP_DR / PROB_SP_DR(Nsp)
end if
call Disp_Rot( MaxSp, Nmol, Nlj, Nion, &
Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, TYPEion, DAMPion, &
BoxSize, DXYZ, DROT, LENGTHlj, LENGTHion, SPECIES, &
STARTlj, STARTion, PROB_DR, PROB_SP_DR, &
Nham, BETA, LNW, LNPSI, LnPi, &
Nljgrs, EPS, SIG, CP, ALP, RMAX, MASSlj, &
Niongrs, CHARGE, MASSion, &
Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, &
ULJBOX, UFOURIER, UREAL, USURF, &
DispOrRot, SpID, Success, Seed )
NATT( DispOrRot, SpID ) = NATT( DispOrRot, SpID ) + 1
if( Success ) NSUC( DispOrRot, SpID ) = NSUC( DispOrRot, SpID ) + 1
MOVETIME( DispOrRot ) = MOVETIME( DispOrRot ) + timef() - tt ! # C_ALPHA
! MOVETIME( DispOrRot ) = MOVETIME( DispOrRot ) + etime(xxx) - tt ! # U_ALPHA
else if( Ranq < PROB_MOVE(2) ) then
tt = timef() ! # C_ALPHA
! tt = etime(xxx) ! # U_ALPHA
if( ran2(Seed) < 0.5 ) then
Ranq2 = ran2(Seed)
SpID = 0
j = 1
do while ( SpID == 0 )
if( Ranq2 < PROB_SP_CD(j) ) SpID = j
j = j + 1
end do
if( tag_val(SpID) == 0 ) then
call Add( Nlj, Nion, MaxMol, MaxSp, MaxBeads, MaxEE, Nmol, &
Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, TYPEion, DAMPion, &
MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
METHOD, MaxInt, MaxReal, INTPARAM, &
REALPARAM, BEADTYPE, BEADDAMP, BONDSAPART, &
ERSTEPS, ERSTART, EREND, EESTEPS, EEW, &
BoxSize, SPECIES, LENGTHlj, LENGTHion, STARTlj, &
STARTion, LENLJ, LENION, GROUPTYPE_LJ, &
GROUPTYPE_ION, Nham, BETA, ZETA, LNW, LNPSI, LnPi, &
XLRCORR, ELRCORR, Nljgrs, EPS, SIG, CP, ALP, RMAX, &
Niongrs, CHARGE, Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, &
ULJBOX, ULJLR, UREAL, UFOURIER, USURF, USELF_CH, &
USELF_MOL, ULJ_MOL, UION_MOL, UINTRA, SpID, &
Success, tag_val(SpID), tag_mol(SpID), &
Nres, Nresmol, reslen, Xresmols, &
Yresmols, Zresmols, Uint_resm, &
Ulj_resm, Uion_resm, &
f14_lj, f14_ion, ER_HIST, Seed )
NATT( 3, SpID ) = NATT( 3, SpID ) + 1
if( Success ) NSUC( 3, SpID ) = NSUC( 3, SpID ) + 1
MOVETIME(3) = MOVETIME(3) + timef() - tt ! # C_ALPHA
! MOVETIME(3) = MOVETIME(3) + etime(xxx) - tt ! # U_ALPHA
else
call alpha_ch( MaxSp, MaxMol, MaxBeads, MaxEE, Nmol, &
Nlj, Nion, Xlj, Ylj, Zlj, &
TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, &
TYPEion, DAMPion, BoxSize, &
LENGTHlj, LENGTHion, SPECIES, &
STARTlj, STARTion, &
BEADTYPE, BEADDAMP, BONDSAPART, &
EESTEPS, EEW, &
Nham, BETA, ZETA, LNW, LNPSI, LnPi, &
Nljgrs, EPS, SIG, CP, ALP, RMAX, &
Niongrs, CHARGE, Alpha, &
Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, &
SUMQX, SUMQY, SUMQZ, ULJBOX, ULJ_MOL, &
UFOURIER, UREAL, USURF, USELF_CH, &
USELF_MOL, UION_MOL, &
tag_val(SpID), tag_mol(SpID), +1, &
f14_lj, f14_ion, &
SpID, Success, Seed )
NATT( 6, SpID ) = NATT( 6, SpID ) + 1
if( Success ) NSUC( 6, SpID ) = NSUC( 6, SpID ) + 1
MOVETIME(6) = MOVETIME(6) + timef() - tt ! # C_ALPHA
! MOVETIME(6) = MOVETIME(6) + etime(xxx) - tt ! # U_ALPHA
end if
else
Ranq2 = ran2(Seed)
SpID = 0
j = 1
do while ( SpID == 0 )
if( Ranq2 < PROB_SP_CD(j) ) SpID = j
j = j + 1
end do
if( tag_val(SpID) == 0 ) tag_val(SpID) = EESTEPS(SpID)
if( tag_val(SpID) == 1 ) then
call Remove( Nlj, Nion, MaxMol, MaxSp, MaxEE, Nmol, &
Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, TYPEion, DAMPion, &
MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, &
REALPARAM, BEADTYPE, BONDSAPART, &
ERSTEPS, ERSTART, EREND, EESTEPS, EEW, &
BoxSize, SPECIES, LENGTHlj, LENGTHion, STARTlj, &
STARTion, Nham, BETA, ZETA, LNW, LNPSI, &
LnPi, XLRCORR, ELRCORR, &
Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, CHARGE, &
Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, &
ULJBOX, ULJLR, UREAL, UFOURIER, USURF, USELF_CH, &
USELF_MOL, ULJ_MOL, UION_MOL, UINTRA, SpID, &
Success, tag_val(SpID), tag_mol(SpID), rmol, &
Nres, Nresmol, reslen, Xresmols, &
Yresmols, Zresmols, Uint_resm, &
Ulj_resm, Uion_resm, &
f14_lj, f14_ion, ER_HIST, Seed )
if( success ) then
do j = 1, Nsp
if( tag_mol(j) > rmol ) then
tag_mol(j) = tag_mol(j) - 1
end if
end do
end if
NATT( 4, SpID ) = NATT( 4, SpID ) + 1
if( Success ) NSUC( 4, SpID ) = NSUC( 4, SpID ) + 1
MOVETIME(4) = MOVETIME(4) + timef() - tt ! # C_ALPHA
! MOVETIME(4) = MOVETIME(4) + etime(xxx) - tt ! # U_ALPHA
else
call alpha_ch( MaxSp, MaxMol, MaxBeads, MaxEE, Nmol, &
Nlj, Nion, Xlj, Ylj, Zlj, &
TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, &
TYPEion, DAMPion, BoxSize, &
LENGTHlj, LENGTHion, SPECIES, &
STARTlj, STARTion, &
BEADTYPE, BEADDAMP, BONDSAPART, &
EESTEPS, EEW, &
Nham, BETA, ZETA, LNW, LNPSI, LnPi, &
Nljgrs, EPS, SIG, CP, ALP, RMAX, &
Niongrs, CHARGE, Alpha, &
Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, &
SUMQX, SUMQY, SUMQZ, ULJBOX, ULJ_MOL, &
UFOURIER, UREAL, USURF, USELF_CH, &
USELF_MOL, UION_MOL, &
tag_val(SpID), tag_mol(SpID), -1, &
f14_lj, f14_ion, &
SpID, Success, Seed )
NATT( 7, SpID ) = NATT( 7, SpID ) + 1
if( Success ) NSUC( 7, SpID ) = NSUC( 7, SpID ) + 1
MOVETIME(7) = MOVETIME(7) + timef() - tt ! # C_ALPHA
! MOVETIME(7) = MOVETIME(7) + etime(xxx) - tt ! # U_ALPHA
end if
end if
else
tt = timef() ! # C_ALPHA
! tt = etime(xxx) ! # U_ALPHA
call ReGrow( Nlj, Nion, MaxMol, MaxSp, Nmol, &
Xlj, Ylj, Zlj, TYPElj, DAMPlj2, DAMPlj3, &
Xion, Yion, Zion, TYPEion, DAMPion, &
NSTEPS, MaxSteps, STEPSTART, STEPLENGTH, NTRIALS, &
MaxBeads, METHOD, MaxInt, MaxReal, INTPARAM, &
REALPARAM, BEADTYPE, BONDSAPART, &
BoxSize, SPECIES, LENGTHlj, LENGTHion, &
STARTlj, STARTion, PROB_SP_RG, &
Nham, BETA, LNW, LNPSI, LnPi, XLRCORR, ELRCORR, &
Nljgrs, EPS, SIG, CP, ALP, RMAX, Niongrs, CHARGE, &
Alpha, Kmax, Nkvec, KX, KY, KZ, CONST, &
EXPX, EXPY, EXPZ, SUMQEXPV, SUMQX, SUMQY, SUMQZ, &
ULJBOX, ULJLR, UREAL, UFOURIER, USURF, USELF_CH, &
USELF_MOL, ULJ_MOL, UION_MOL, UINTRA, SpID, Success, &
Nres, Nresmol, reslen, Xresmols, &
Yresmols, Zresmols, Uint_resm, &
Ulj_resm, Uion_resm, &
f14_lj, f14_ion, Seed )
NATT( 5, SpID ) = NATT( 5, SpID ) + 1
if( Success ) NSUC( 5, SpID ) = NSUC( 5, SpID ) + 1
MOVETIME(5) = MOVETIME(5) + timef() - tt ! # C_ALPHA
! MOVETIME(5) = MOVETIME(5) + etime(xxx) - tt ! # U_ALPHA
end if
do h = 1, Nham
UION(h) = USURF(h) + UFOURIER(h) + UREAL(h) - USELF_CH(h) - &
sum( USELF_MOL( 1:Nmol(0), h) ) + &
sum( UION_MOL( 1:Nmol(0), h) )
ULJ(h) = ULJBOX(h) + ULJLR(h) + &
sum( ULJ_MOL( 1:Nmol(0), h) )
UTOT(h) = UION(h) + ULJ(h)
end do
do j = 1, Nsp
if( tag_val(j) == 0 ) then
if( Nmol(j) > nt1(j) ) nin(j) = nin(j) + 1
if( Nmol(j) < nt1(j) ) nout(j) = nout(j) + 1
nt1(j) = Nmol(j)
end if
end do
! for EE weights
do k = 1, Nsp
ALP_PROB(tag_val(k),k) = ALP_PROB(tag_val(k),k) + 1
end do
! for HS stuff
PSI_PI = exp( LNPSI - LnPi )
AVG_PSI_PI = AVG_PSI_PI + PSI_PI
if( i > StartConf + sum( NST_TOT(1:stage) ) - NST_TOT(stage) + NST_WT(stage) ) then ! For calculation of W
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -