📄 gcmc.f90
字号:
U_AVG = U_AVG + UTOT
PSI_PI_AVG = PSI_PI_AVG + PSI_PI
WCounts = WCounts + 1
end if
if( .NOT. FromDisk .AND. i > int( 0.25 * Nequil ) &
.AND. i <= Nequil ) then
do j = 1, Nsp
Nsim_min(j) = min( Nmol(j), Nsim_min(j) )
Nsim_max(j) = max( Nmol(j), Nsim_max(j) )
end do
Usim_min = min( minval(UTOT), Usim_min )
Usim_max = max( maxval(UTOT), Usim_max )
else if( FromDisk .AND. i <= StartConf + Nequil ) then
do j = 1, Nsp
Nsim_min(j) = min( Nmol(j), Nsim_min(j) )
Nsim_max(j) = max( Nmol(j), Nsim_max(j) )
end do
Usim_min = min( minval(UTOT), Usim_min )
Usim_max = max( maxval(UTOT), Usim_max )
end if
if( mod( i, NinCycle ) == 0 ) then
write(*,*)
write(*,"(A,T15,I8)") ' MC Step = ', i
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)
! format 11 is format( A, T15, <Nham>F15.4 )
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)
write(*,"(A,T15,4I8)") ' n in = ', nin(1:Nsp)
write(*,"(A,T15,4I8)") ' n out = ', nout(1:Nsp)
write(*,"(A,T15,10I8)") ' probs = ', ALP_PROB(0:EESTEPS(1)-1,1:Nsp)
tot_time = ( timef() - tstart ) ! # C_ALPHA
! tot_time = ( etime(xxx) - tstart ) ! # U_ALPHA
if ( tot_time < 60.0 ) then
write(*, "(1X, A, T15, F15.4, A)") &
'Total Time = ', tot_time, ' sec.'
else if ( tot_time < 3600.0 ) then
write(*, "(1X, A, T15, F15.4, A)") &
'Total Time = ', tot_time/60.0, ' min.'
else
write(*, "(1X, A, T15, F15.4, A)") &
'Total Time = ', tot_time/3600.0, ' hrs.'
endif
open( 50, file = PlotsFile, position = 'append' )
write(50,12) i, &
dot_product( Nmol(1:Nsp) , MASSsp(1:Nsp) ) * 1.0e27 / Nav / ( Volume ), &
UTOT(:), LNPSI(:), Nmol(1:Nsp)
12 format( I8, <2*Nham+1>(1X,G15.7), <Nsp>I8 )
close(50)
Moves = 0
end if
if( mod( i, Nrefresh ) == 0 ) then
do m = 1, Nres
do j = 1, Nresmol
Xres(1:reslen(m),m) = Xresmols(j,1:reslen(m),m)
Yres(1:reslen(m),m) = Yresmols(j,1:reslen(m),m)
Zres(1:reslen(m),m) = Zresmols(j,1:reslen(m),m)
do k = 1, Nrefrmoves
call ResUpdate( reslen(m), Xres(:,m), &
Yres(:,m), Zres(:,m), &
Nresmv(m), MaxResmv, MaxResp, &
PROB_RMV(:,m), RESMOVE(:,m), &
RESPARAM(:,:,m), &
Nvib(m), Nbend(m), Ntor(m), &
MaxVib, MaxBend, MaxTor, &
IVIB(:,:,m), IBEND(:,:,m), ITOR(:,:,m), &
RVIB(:,:,m), RBEND(:,:,m), RTOR(:,:,m), &
RES_ATT, RES_SUC, &
MaxBeads, &
BEADTYPE(:,m), GROUPTYPE(:,m), &
BONDSAPART(:,:,m), Nham, Nljgrs, &
EPS, SIG, CP, ALP, RMAX, &
Niongrs, CHARGE, &
f14_lj, f14_ion, &
BETA, LNW, BoxSize, &
UResInt(m), UResLj(:,m), UResIon(:,m), &
Seed, th_h, phi_h )
end do
Xresmols(j,1:reslen(m),m) = Xres(1:reslen(m),m)
Yresmols(j,1:reslen(m),m) = Yres(1:reslen(m),m)
Zresmols(j,1:reslen(m),m) = Zres(1:reslen(m),m)
Uint_resm(j,m) = UResInt(m)
Ulj_resm(j,:,m) = UResLj(:,m)
Uion_resm(j,:,m) = UResIon(:,m)
end do
end do
end if
if( mod( i, Ncollect ) == 0 .AND. i > StartConf + Nequil &
.AND. sum(tag_val) == 0 ) then
PSI_PI = PSI_PI * exp( LNW ) * real( Nham )
do h = 1, Nham
do j = 1, Nsp
N_HIST( Nmol(j), j, h ) = N_HIST( Nmol(j), j, h ) + PSI_PI(h)
Nsim_min(j) = min( Nmol(j), Nsim_min(j) )
Nsim_max(j) = max( Nmol(j), Nsim_max(j) )
end do
end do
if( histtype == 'list' ) then
nentry = nentry + 1
do j = 1, Nsp
Npart( nentry, j ) = Nmol(j)
end do
Energy( nentry ) = UTOT(1)
Usim_min = min( UTOT(1), Usim_min )
Usim_max = max( UTOT(1), Usim_max )
else if( histtype == 'matrix' ) then
do h = 1, Nham
if( UTOT(h) > Umin .AND. UTOT(h) < Umax ) then
j = int( ( UTOT(h) - Umin ) / Uwidth ) + 1
UN_HIST( j, Nmol(0), h ) = UN_HIST( j, Nmol(0), h ) + PSI_PI(h)
else
if( UTOT(h) < Umin ) then
write(*,*) ' Umin was exceeded! '
write(*,*) ' Umin = ', Umin
write(*,*) ' Energy = ', UTOT(h)
else if( UTOT(h) > Umax ) then
write(*,*) ' Umax was exceeded! '
write(*,*) ' Umax = ', Umax
write(*,*) ' Energy = ', UTOT(h)
end if
write(*,*) ' The program will stop. '
stop
end if
Usim_min = min( UTOT(h), Usim_min )
Usim_max = max( UTOT(h), Usim_max )
end do
end if
end if
if( mod( i, Nadjust ) == 0 .AND. i <= StartConf + Nequil ) then
TOT_NATT = TOT_NATT + NATT
TOT_NSUC = TOT_NSUC + NSUC
call NewMaxima( Nsp, MaxSp, nmtype, DXYZ, DROT, NATT, NSUC )
end if
if( mod( i, Neew ) == 0 .AND. i <= StartConf + Nequil ) then
do j = 1, Nsp
call EEWeights( EESTEPS(j), EEW(:,j), ALP_PROB(0:EESTEPS(j),j) )
write(*,*)
write(*,"(A,I2)") ' New EE Weights for Species ', j
write(*,"(10f10.6)") EEW(1:EESTEPS(j),j)
write(*,*)
end do
ALP_PROB = 0
end if
if( mod( i, Nstorage ) == 0 ) then
if( i <= StartConf + Nequil ) then
call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
i, DXYZ, DROT, MaxMol, Ubins, Usim_min, Usim_max, &
tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
Nsim_min(1), Nsim_max(1), Nham, LNW, Seed )
else
if( histtype == 'matrix' ) then
call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
i, DXYZ, DROT, MaxMol, Ubins, Umin, Umax, &
tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
Nsim_min(1), Nsim_max(1), Nham, LNW, Seed )
call WriteHist( MaxMol, MaxSp, Ubins, Nham, BETA, ZETA, Uwidth, Umin, &
UN_HIST, N_HIST, BoxSize, UN_HistFile, N_HistFile(1), &
Nsim_min, Nsim_max )
else if( histtype == 'list' ) then
call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
i, DXYZ, DROT, MaxMol, Ubins, Usim_min, Usim_max, &
tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
Nsim_min(1), Nsim_max(1), Nham, LNW, Seed )
open(90, file=UN_HistFile, position='append')
do j = 1, nentry
write(90,14) Npart(j,1:Nsp), Energy(j)
14 format(1x,<Nsp>I4,G16.7)
end do
close(90)
nentry = 0
do j = 1, MaxSp
open(100, file=N_HistFile(j))
do k = Nsim_min(j), Nsim_max(j)
write(100,"(I8,G16.7E3)") k, N_HIST(k,j,1)
end do
close(100)
end do
end if
end if
open(19, file = 'er.dat')
do k = 1, Nsp
do j = 1, ERSTEPS(1) + 1
write(19,*) j, ER_HIST(j, k, 1:2)
end do
write(19,*)
end do
close(19)
end if
if( i == StartConf + sum( NST_TOT(1:stage) ) ) then
open(90, file = WFile, position = 'append')
write(90,*)
write(90,"(A,I4)") ' Results after stage ', stage
write(90,*)
if( STOP_H(stage) < 0 ) then
STOP_H(stage) = - STOP_H(stage)
! U_AVG = U_AVG / real( WCounts )
! do h = 2, Nham
! BETA(h) = BETA(1) * U_AVG(1) / U_AVG(h)
! end do
call predict( Nsim_min(1), Nsim_max(1), Ubins, Umin, Uwidth, Nham, &
MaxMol, UN_HIST, BETA, ZETA )
if( STOP_H(stage) /= 100 ) then
PSI_PI_AVG = PSI_PI_AVG / real( WCounts )
call Weights( Nham, STOP_H(stage), LNW, PSI_PI_AVG )
end if
write(90,"(A)") ' Hamiltonian Weight Temperature Activity '
do h = 1, Nham
write(90,"(5x,I3,T15,G14.7E3,T32,F9.4,T46,F8.4)") h, exp( LNW(h) ), 1.0/BETA(h), log(ZETA(1,h))
end do
else if( STOP_H(stage) > 0 ) then
LNW_old = LNW
PSI_PI_AVG = PSI_PI_AVG / real( WCounts )
call Weights( Nham, STOP_H(stage), LNW, PSI_PI_AVG )
write(90,"(A)") ' Hamiltonian Weight PSI_PI_AVG Coverage '
do h = 1, Nham
write(90,"(5x,I3,T15,G14.7E3,T32,G14.7E3,T47,F8.5)") &
h, exp( LNW(h) ), PSI_PI_AVG(h), PSI_PI_AVG(h) * exp(LNW_old(h)) * real(Nham)
end do
end if
close(90)
PSI_PI_AVG = 0.0
WCounts = 0
if( i > StartConf + Nequil ) then
UN_HIST = 0.0
N_HIST = 0.0
end if
Nsim_min = MaxMol
Nsim_max = -1
LNPSI(:) = Nmol(0)*log(ZETA(1,:)) + 3.0*Nmol(0)*log(BoxSize) - BETA(:)*UTOT(:)
do j = 1, Nsp
do k = 1, Nmol(j)
LNPSI = LNPSI - log(real(k))
end do
end do
Largest = maxval( LNW + LNPSI )
LnPi = log( sum( exp( LNW + LNPSI - Largest ) ) ) + Largest
stage = stage + 1
if( Stage > Nstages ) stage = Nstages
end if
if( i == startConf + Nequil .OR. i == StartConf + Nequil + Nprod ) then
TOT_NATT = TOT_NATT + NATT
TOT_NSUC = TOT_NSUC + NSUC
if( i == StartConf + Nequil ) then
EquilOrProd = 1
stc = StartConf + 1
else
EquilOrProd = 2
stc = StartConf + Nequil + 1
end if
if( EquilOrProd == 1 ) then
if( histtype == 'matrix' ) then
Urange = Usim_max - Usim_min
Umax = Usim_max + SafetyFactor * Urange
Umin = Usim_min - SafetyFactor * Urange
if( Umax > 0 ) Umax = ( real( int( Umax / Uwidth ) ) + 1.5 ) * Uwidth
if( Umax < 0 ) Umax = ( real( int( Umax / Uwidth ) ) + 0.5 ) * Uwidth
if( Umin > 0 ) Umin = ( real( int( Umin / Uwidth ) ) - 0.5 ) * Uwidth
if( Umin < 0 ) Umin = ( real( int( Umin / Uwidth ) ) - 1.5 ) * Uwidth
Ubins = nint( ( Umax - Umin ) / Uwidth )
allocate( UN_HIST( Ubins, 0:MaxMol, Nham ) )
UN_HIST = 0.0
end if
end if
tot_time = timef() - tstart ! # C_ALPHA
! tot_time = etime(xxx) - tstart ! # U_ALPHA
call WriteResults( ResultsFile, EquilOrProd, stc, i, &
Nsp, MaxSp, nmtype, TOT_NATT, TOT_NSUC, &
NAMEsp, MOVETIME, tot_time, &
MaxEE, EESTEPS, ALP_PROB, nin, nout, &
Usim_min, Usim_max, &
Umin, Umax, Ubins, Nsim_min, Nsim_max )
if( EquilOrProd == 1 .AND. Nprod > 0 ) then
Nsim_min = MaxMol
Nsim_max = 0
Usim_min = Umax
Usim_max = Umin
end if
moves = 0
NATT = 0
NSUC = 0
TOT_NATT = 0
TOT_NSUC = 0
ALP_PROB = 0
nin = 0
nout = 0
MOVETIME = 0.0
end if
end do
Istore = 0
i = StartConf + Nequil + Nprod
call StoreConf( Istore, FinalConf, Nsp, MaxSp, Nmol, NAMEsp, &
i, DXYZ, DROT, MaxMol, Ubins, Umin, Umax, &
tag_val, tag_mol, MaxEE, EESTEPS, EEW, &
LENGTHlj, LENGTHion, SPECIES, MaxBeads, BEADTYPE, &
MaxLJ, MaxIon, Xlj, Ylj, Zlj, Xion, Yion, Zion, &
Nsim_min, Nsim_max, Nham, LNW, Seed )
write(*,*)
write(*,11) ' UTotal = ', UTOT(1:Nham)
write(*,11) ' LNPSI = ', LNPSI(1:Nham)
write(*,*)
close(50)
end program gcmc
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -