|
|
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
zh1987hs(金币+10):谢谢 2010-12-04 23:13:56
ghcacj(金币+3, 模拟EPI+1):精彩答疑 2010-12-06 12:16:56
答应过的事情, 虽然觉得其实对楼主帮助也不大. 因为我不可能给你写一个完整的GibbsMC, 如果你自己不努力去理解, 那么什么都是白搭.
这是大体结构. 按照需要填充就是了. 应该此贴不会再更新了.
文件名: GibbsMC.f90
文件内容如下:
PROGRAM MC
USE global
IMPLICIT NONE
! Local
INTEGER :: ibox, iblock, istep
INTEGER :: ts,tf,rate
DOUBLE PRECISION :: engtst, virtst, engerr
! Initialize the random number generator
CALL seed_random(idum)
! Read the input file
CALL read_input
! Get the start time
CALL system_clock(count=ts, count_rate=rate)
! Perform minimization
CALL relax
! Calculate energy and virial
DO ibox = 1,n_box
CALL eng_total(ibox,energy(ibox),virial(ibox))
CALL write_config(ibox,0)
ENDDO
! Loop over the blocks
DO iblock = 1,n_blocks
! Reset basic statistics
CALL stats(1,0,0)
! Reset optional statistics
IF(opt .EQ. 1) CALL sample(1,0,0)
! Loop over the steps per block
DO istep = 1,n_steps
! Perform a trial move
CALL trial_move(iblock,istep)
! Accumulate basic statistics
CALL stats(2,iblock,istep)
! Accumulate optional statistics
IF(opt .EQ. 1) call sample(2,iblock,istep)
! End loop over the steps per block
ENDDO
! Average basic statistics for iblock
CALL stats(3,iblock,0)
! Average optional statistics for iblock
IF(opt .EQ. 1) call sample(3,iblock,0)
! Update the max values (disp, volchng, etc)
CALL update_max
! Print status
IF(iblock .LE. n_equil) THEN
IF(iblock .EQ. 1) WRITE(*,'(A)')'Equilibration stage... '
IF(MOD(NINT(100.0d0*DBLE(iblock)),10*n_equil) .EQ. 0)&
&WRITE(*,'(I0,A3)')NINT(100.0d0*DBLE(iblock)/DBLE(n_equil)),'% '
IF(iblock .EQ. n_equil) WRITE(*,*)''
ELSE
IF(iblock .EQ. n_equil + 1) WRITE(*,'(A)')'Production stage... '
IF(MOD(NINT(100.0d0*DBLE(iblock-n_equil)),10*(n_blocks-n_equil)).EQ. 0)&
&WRITE(*,'(I0,A3)')NINT(100.0d0*DBLE(iblock-n_equil)/DBLE(n_blocks-n_equil)),'% '
IF(iblock .EQ. n_blocks) WRITE(*,*)''
ENDIF
! Check the total and updated energy consistency for each box (every 10 blocks)
IF(MOD(iblock,1) .EQ. 0) THEN
DO ibox = 1,n_box
CALL eng_total(ibox,engtst,virtst)
engerr = dABS((engtst - energy(ibox)))
IF (engerr .GE. 1.0d-5) THEN
WRITE(*,*) 'WARNING: LARGE ERROR BETWEEN THE TOTAL AND UPDATED ENERGY'
WRITE(*,'(A7,I0,A6,I0,A10,F20.8,A10,F20.8)')&
&'BLOCK: ',iblock,' BOX: ',ibox,' UPD ENG: ',energy(ibox),' TOT ENG: ',engtst
ENDIF
ENDDO
ENDIF
! End loop over the blocks
ENDDO
! Get the finish time
CALL system_clock(count=tf, count_rate=rate)
! Write the elapsed time to screen
CALL time_elapsed(tf,ts,rate)
! Write basic statistics
CALL stats(4,0,0)
! Write optional statistics
IF(opt .EQ. 1) CALL sample(4,n_blocks+1,0)
! Write the final configurations
DO ibox = 1,n_box
CALL write_config(ibox,n_blocks+1)
ENDDO
END PROGRAM MC |
|
回复此楼
点击这里搜索更多相关资源
20