应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Monte Carlo » 【原创】抛土引玉, Monte Carlo 如何自己写程序: NVT_LJ
811/9返回列表123456下一页末页
查看: 8713  |  回复: 80
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 2 个 模拟EPI ，点击这里进行查看

qphll

金虫 (正式写手)

[交流] 【原创】抛土引玉, Monte Carlo 如何自己写程序: NVT_LJ已有57人参与

引子: 注册小木虫很多年了, 一直是条小虫子, 也喜欢这样的论坛, 有虫, 有大虫, 有菜, 也有肉. 在学术和科研这个大盘子上面, 谁都是小虫子, 谁都能成为某一方向的专家.或许大家内心偶尔还是有那种"为中华之崛起而读书"的热血冲动, 但是更多时候, 为现实的残酷而黯然, 为学海之无涯而沮丧. 但是正像虫卵可以慢慢长成小虫一样, 只要你不断努力, 小虫也会破茧成为美丽的蝴蝶.

不管哭与笑, 不论败与成, Keep fighting, Never give up.

与大家共勉.

qphll
Nov. 24, 2010
回复此楼

» 收录本帖的淘帖专辑推荐

材料计算模拟实用技巧 分子模拟 生活休闲 学习方法
condensed matter physics 第一性原理和电化学 量子化学理论及量化软件专辑 量化,第一性,MD笔记
分子生物物理理论

» 本帖已获得的红花（最新10朵）

honglyn zhang_jaj nilelixp

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
Life, Love, Laugh.
1楼2010-11-25 00:40:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):鼓励 2010-11-25 09:23:58
言归正传, 这一个帖子的目的, 是分享与讨论. 大家习惯于Material Studio的方便, 也赞赏GNU下各个package的免费分享. 但是严肃的计算是建立在透彻理解原理的基础上的, 这也是为什么就算限于规模和条件, 只能使用别人的软件, 你至少也应该看看源代码中关键的部分.

关于Monte Carlo, 不用多说, 大家知道那是一个好东西. 要不然为啥美帝两弹的某部分计算需要依靠它呢? 传统意义上的分子模拟, 也是从五十年代的第一篇Monte Carlo文章拉开序幕, 闪亮登场的.

但是大家, 特别是作为我们学生, 不是那么喜欢Monte Carlo, 因为.... 要自己写程序.
汪老师翻译的书, 多少年都是宝典. 确实在Monte Carlo那块, 也没有什么书能超越那本. 但是大家想过没有, 这其实也是暗示着, Monte Carlo那边, 核心的东西也就是那么一些, 颇有些"一旦拥有, 万世无愁"的意味. 大家数数吧, 有多少人, 多少课题组的研究, 是依靠一些也许在十年前写的Monte Carlo代码? 当然算法的发展是有的, 可是, 你明白的, 一旦有了核心的代码, 修修补补是容易的.

我这里贴出来一个完整的Monte Carlo代码, 是关于NVT下面Argon的Monte Carlo.

每行代码都是简单易懂的. 当然, 我也没有做多少优化, 几乎是采用白话文的方式写下来练手. 对你我而言, 建立起自己可以写代码的信心, 是压倒性地头等重要. 从这个源码, 你知道, "我能行!!!"

好, 上菜!!
赞一下(6人) 回复此楼
Life, Love, Laugh.
3楼2010-11-25 00:56:36
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10, 模拟EPI+1):鼓励 2010-11-25 09:24:33
MODULE global_variables
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Contains all the global variables that will be used in the program !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  !     qphll                                                          !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Number of particles  
!----------------------------------------------------------------------!  
  ! Number of nitrogen molecules
  INTEGER :: n_molecules
  
  INTEGER,PARAMETER :: n_molecules_max=1000
  

! Positions
!----------------------------------------------------------------------!
  ! Positions of the centers of mass of the N molec.
  DOUBLE PRECISION,DIMENSION(n_molecules_max) :: rx, ry, rz
  

! Geometry  
!----------------------------------------------------------------------!
  ! Simulation box length
  DOUBLE PRECISION :: box

! Lennard Jones parameters
  DOUBLE PRECISION :: sig
  DOUBLE PRECISION :: eps

! Physical variables
!----------------------------------------------------------------------!
  ! Temperature (non-dimensional)
  DOUBLE PRECISION :: temp

  ! Energy, etc
  DOUBLE PRECISION :: energy

  ! Block accumulators
  DOUBLE PRECISION :: blkener, blkpress, blkenersq

  ! Virial and pressure
  DOUBLE PRECISION :: pressure

! Parameters
!----------------------------------------------------------------------!
  ! Cutoff distances
  DOUBLE PRECISION :: cutoff_nn, cutoff_nn_square

  ! Very small distances will cause huge energies
  DOUBLE PRECISION, PARAMETER :: huge_number = 10000
  DOUBLE PRECISION :: small_rsq_nn, huge_energy


! Acceptance/Rejections Statistics
!----------------------------------------------------------------------!
  ! Translation statistics
  INTEGER :: acctrans, ntrans


! Monte Carlo parameters
!----------------------------------------------------------------------!
  ! Number of blocks
  INTEGER :: nblocks

  ! Number of steps per block
  INTEGER :: nsteps

  ! Maximum displacement in a MC move
  DOUBLE PRECISION :: maxdisp

  ! Adjust the parameter maxdisp every nadjust blocks
  INTEGER :: nadjust

! Radial distribution function
!----------------------------------------------------------------------!
  INTEGER,PARAMETER :: maxbin = 200
  INTEGER, DIMENSION(maxbin) :: hist
  DOUBLE PRECISION :: delr,constant
  INTEGER :: nsamplegr,totalsamplesgr

! DCD file output
!----------------------------------------------------------------------!
  INTEGER :: dcd_freq

! Others
!----------------------------------------------------------------------!
  ! Dummy variable for random number generator
  INTEGER :: idum

  ! Random number generator
  DOUBLE PRECISION,EXTERNAL :: random

  ! Pi
  DOUBLE PRECISION,PARAMETER :: Pi = 3.141592d0
  DOUBLE PRECISION,PARAMETER :: Kb = 1.3806d-23  ! J/K
  DOUBLE PRECISION,PARAMETER :: Na = 6.0221d23   ! 1/mol
  DOUBLE PRECISION,PARAMETER :: R = Kb*Na        ! J/(mol K)
  DOUBLE PRECISION,PARAMETER :: Pp = Kb*1.0d+30  ! (J cu.A)/(K cu.M) OR (Pa cu.A)/K
  DOUBLE PRECISION,PARAMETER :: Pb = Pp*1.0d-5   ! (bar cu. A)


END MODULE global_variables

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM nvt

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Main program                                                       !
  ! Procedures used:                                                   !
  !  - start, runner, finish                                           !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  CALL start
  CALL runner
  CALL finish

  STOP

END PROGRAM nvt

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE runner

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Performs 'nblocks' blocks of 'nsteps' Monte Carlo steps            !
  !                                                                    !
  ! Procedures used:                                                   !
  !  - random                                                          !
  !  - translate                                                       !
  !  - total_energy                                                    !
  !  - adjust_mc_parameters                                            !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Local variables
  INTEGER :: iblock, istep
  DOUBLE PRECISION :: etest, ptest, Utest, Wtest

  ! Calculate the total energy
  CALL total_energy(energy,pressure)

  ! Perform runs over blocks
  DO iblock=1,nblocks

     ! Set block accumulators to zero
     blkener = 0.0d0
     blkenersq = 0.0d0
     blkpress = 0.0d0

     ! Perform runs over steps
     DO istep=1,nsteps

        ! Perform an MC step
        CALL translate

        ! Update block accumulators
        blkener = blkener + energy
        blkpress = blkpress + pressure
        blkenersq = blkenersq + energy**2

        ! End loop over steps
     END DO


        ! Average properties at the end of the block
        blkener = blkener / DBLE(nsteps)/DBLE(n_molecules)
        blkenersq = blkenersq / REAL(nsteps)/DBLE(n_molecules)**2
        blkpress = blkpress / DBLE(nsteps)

    CALL total_energy(energy,pressure)

        ! Write bulk properties to file
        IF(sig == 1.0 .AND. eps == 1.0) THEN

          ! Write results in output file
          WRITE(40,FMT='(I5,F18.6,F24.6,F18.6)') &
               &iblock,blkener,blkenersq,&
               &((n_molecules*temp-blkpress)/box**3)

        ELSE

          ! Write results in output file
          WRITE(40,FMT='(I5,F18.6,F24.6,F18.6)') &
               &iblock,R*blkener,(R**2)*blkenersq,&
               &Pb*((n_molecules*temp-blkpress)/box**3)

        END IF


        ! Adjust MC parameters
        IF (MOD(iblock,nadjust) == 0) CALL adjust_mc_parameters

        ! Sample radial distribution function
        IF (MOD(iblock,nsamplegr) == 0) CALL sample_gr

        ! Write to the DCD file
        CALL write_dcd(iblock)
   
     ! End loop over blocks
  END DO

  ! Close output file
  CLOSE(UNIT=40)

  RETURN

END SUBROUTINE runner

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



REAL FUNCTION random(idum)

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Generates a random number between 0.0 and 1.0                      !
  ! Modification of ran3, taken from "Numerical Recipes in Fortran"    !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  IMPLICIT NONE

  
  ! Local variables
  REAL, PARAMETER :: mbig=4000000., mseed=1618033., mz=0., fac=1./mbig
  INTEGER :: idum,i, ii,  k
  INTEGER, SAVE :: iff, inext, inextp
  REAL :: mj, mk
  REAL, SAVE :: ma(55)
  DATA iff /0/


  ! Initialization
  IF(idum<0.OR.iff==0) THEN
     iff=1
     mj=mseed-IABS(idum)
     mj=MOD(mj, mbig)
     ma(55)=mj
     mk=1

     DO i=1,54
        ii=MOD(21*i,55)
        ma(ii)=mk
        mk=mj-mk
        IF(mk         mj=ma(ii)
     END DO
     DO k=1,4
        DO i=1,55
           ma(i)=ma(i)-ma(1+MOD(i+30,55))
           IF(ma(i)         END DO
     END DO
     inext=0
     inextp=31
     idum=1
  END IF

  ! Start to generate the random number

  inext=inext+1
  IF(inext==56) inext=1
  inextp=inextp+1
  IF(inextp==56) inextp=1
  mj=ma(inext)-ma(inextp)
  IF(mj   ma(inext)=mj
  random=mj*fac
  RETURN

END FUNCTION random

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE  initialize_idum(idum)

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Initializes the seed of the random number generator                !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  IMPLICIT NONE

  ! Local variables
  INTEGER :: idum
  INTEGER, DIMENSION(8) :: val

  idum = -478

!!$  CALL DATE_AND_TIME(VALUES=val)
!!$  idum=-(1000*val(7)+val(8))
  RETURN

END SUBROUTINE initialize_idum

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


SUBROUTINE translate

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Generates a random displacement of a molecule itest                !
  ! selected randomly                                                  !
  !                                                                    !
  ! Procedures used:                                                   !
  !  - random(idum): random number generator                           !
  !  - energy_calculation                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Local variables
  INTEGER :: imol
  DOUBLE PRECISION :: rx_old, ry_old, rz_old
  DOUBLE PRECISION :: rx_new, ry_new, rz_new
  DOUBLE PRECISION :: vold, vnew, deltav,virialold,virialnew
  LOGICAL :: accepted


  ! Select the particle to be translated
  imol = INT(random(idum) * n_molecules) + 1
  
  ! Read the current position of the molecule itest
  rx_old = rx(imol)
  ry_old = ry(imol)  
  rz_old = rz(imol)

  ! Calculate energy with the old position
  CALL energy_calculation(imol,rx_old,ry_old,rz_old,vold,virialold)

  ! Assign the new position randomly
  ! The maximum displacement is given by maxdisp
  rx_new = rx_old + maxdisp * (2.0d0*random(idum)-1.0d0)
  ry_new = ry_old + maxdisp * (2.0d0*random(idum)-1.0d0)
  rz_new = rz_old + maxdisp * (2.0d0*random(idum)-1.0d0)

  ! Apply periodic boundary conditions
  rx_new = rx_new - DNINT(rx_new/box)*box
  ry_new = ry_new - DNINT(ry_new/box)*box
  rz_new = rz_new - DNINT(rz_new/box)*box

  ! Calculate the energy with the new position
  CALL energy_calculation(imol,rx_new,ry_new,rz_new,vnew,virialnew)

  ! Calculate difference in potential energy
  deltav = vnew - vold

  ! Apply the acceptance criteria
  IF (deltav < 0.0d0) THEN
     accepted = .TRUE.
  ELSEIF (deltav >= (huge_energy-1.0d0)) THEN
     accepted = .FALSE.
  ELSEIF (random(idum) < dEXP(-deltav/temp)) THEN
     accepted = .TRUE.
  ELSE
     accepted = .FALSE.
  END IF

  ! If the move is accepted...
  IF (accepted) THEN
     
     ! Update the number of accepted translation moves
     acctrans = acctrans + 1

     ! Update positions
     rx(imol) = rx_new
     ry(imol) = ry_new
     rz(imol) = rz_new

     ! Update the value of the energy
     energy = energy + deltav

     ! Update the value of pressure
     pressure = pressure + virialnew - virialold
  END IF

  ! Update the number of translation moves
  ntrans = ntrans + 1
!WRITE(*,10) imol, vold, vnew, energy, pressure
!10 FORMAT(I3, 2X F11.6, 2X, F11.6, 2X, F11.6, 2X, F11.6, 2X, F11.6)

  RETURN

END SUBROUTINE translate
赞一下(2人) 回复此楼
Life, Love, Laugh.
4楼2010-11-25 00:58:35
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

香薷谷

金虫 (小有名气)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
虽然 不是做蒙特卡洛，但是为楼主的一句“为中华崛起而读书”感动~顶一下！
赞一下(2人) 回复此楼
43楼2011-01-11 23:05:54
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

zyj8119

木虫 (著名写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
引用回帖:
Originally posted by qphll at 2010-11-25 00:40:07:
引子: 注册小木虫很多年了, 一直是条小虫子, 也喜欢这样的论坛, 有虫, 有大虫, 有菜, 也有肉. 在学术和科研这个大盘子上面, 谁都是小虫子, 谁都能成为某一方向的专家.或许大家内心偶尔还是有那种"为中华之崛 ...

NVT-MC貌似只能看FRENKEL的那本书了，那本书对于系综讲解的很透彻。
赞一下(1人) 回复此楼
好好学习，天天向上。
2楼2010-11-25 00:46:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
(续前, 没办法, 有长度限制....)


SUBROUTINE energy_calculation(itest,rxi,ryi,rzi,potenergy,virial)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Calculates the potential energy of molecule"i" (test molecule)     !
  ! with all the other molecules/atoms ("j"                           !
  !                                                                    !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Variables with argument association
  INTEGER :: itest
  DOUBLE PRECISION :: rxi,ryi,rzi,potenergy,virial


  ! Local variables
  DOUBLE PRECISION :: energy_ij,virial_ij
  DOUBLE PRECISION :: factor
  DOUBLE PRECISION :: rxij,ryij,rzij,rijsq
  INTEGER :: j


  ! Contribution due to interactions with other nitrogen molecules


  !Initialize variables
  potenergy = 0.0d0
  virial = 0.0d0

  ! Start loop over all nitrogen molecules
  DO j=1,n_molecules

     ! avoid self-interactions
     IF (j==itest) CYCLE

     ! calculate vector between centers of mass
     rxij = rx(j) - rxi
     ryij = ry(j) - ryi
     rzij = rz(j) - rzi

     ! apply minimum image convention
     rxij = rxij - DNINT(rxij/box)*box
     ryij = ryij - DNINT(ryij/box)*box
     rzij = rzij - DNINT(rzij/box)*box

     ! Calculate the distance between centers of mass
     rijsq = rxij**2 + ryij**2 + rzij**2

     ! Apply cutoff distance
     IF (rijsq <= cutoff_nn_square) THEN

        ! Calculate the potential due to dispersion interactions

        ! Check if the distance is very small
        IF (rijsq <= small_rsq_nn) THEN
           potenergy = huge_energy
           RETURN
        END IF

        ! Evaluate the L-J potential
        factor = (sig**2/ rijsq)**3
        energy_ij = factor*(factor-1.0d0)

        ! Evaluate the virial
        virial_ij = factor*(1.0d0-2.0d0*factor)

        ! Add contributions to energy and virial
        potenergy = potenergy + energy_ij
        virial = virial + virial_ij

        ! End the calculations for distance smaller than cutoff
     END IF

     ! End loop over nitrogen molecules
  END DO

  ! Multiply energy and virial by the L-J factor
  potenergy = 4.0d0 * eps * potenergy
  virial = 8.0d0 * eps * virial


END SUBROUTINE energy_calculation

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


SUBROUTINE total_energy(potenergy,virial)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Calculates the total energy of the system                          !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Variables with argument association
  DOUBLE PRECISION :: potenergy,virial

  ! Local variables
  DOUBLE PRECISION :: rxi,ryi,rzi
  DOUBLE PRECISION :: energy_ij, virial_ij
  DOUBLE PRECISION :: factor
  DOUBLE PRECISION :: rxij,ryij,rzij,rijsq
  INTEGER :: j,imol


  ! Initialize parameters
  potenergy = 0.0d0
  virial = 0.0d0

  ! Start loop over all nitrogen molecules
  DO imol=1,n_molecules-1

     rxi = rx(imol)
     ryi = ry(imol)
     rzi = rz(imol)


     ! Contribution due to interactions with other nitrogen molecules


     ! Start loop over other nitrogen molecules
        DO j=imol+1,n_molecules

           ! calculate vector between centers of mass
           rxij = rx(j) - rxi
           ryij = ry(j) - ryi
           rzij = rz(j) - rzi

           ! apply minimum image convention
           rxij = rxij - DNINT(rxij/box)*box
           ryij = ryij - DNINT(ryij/box)*box
           rzij = rzij - DNINT(rzij/box)*box

           ! Calculate the distance between centers of mass
           rijsq = rxij**2 + ryij**2 + rzij**2

           ! Apply cutoff distance
           IF (rijsq <= cutoff_nn_square) THEN

              ! Calculate the potential due to dispersion interactions

              ! Check if the distance is very small
              IF (rijsq <= small_rsq_nn) THEN
                 energy = huge_energy
                 PRINT*, 'ERROR: Try with a different initial config'
                 STOP
              END IF

              ! Evaluate the L-J potential
              factor = (sig**2/ rijsq)**3
              energy_ij = factor*(factor-1.0d0)

              ! Calculate the virial
              virial_ij = factor*(1.0d0-2.0d0*factor)

              ! Add contributions to energy and pressure
              potenergy = potenergy + energy_ij
              virial = virial + virial_ij

              ! End the calculations for distance smaller than cutoff
           END IF

           ! End loop over nitrogen molecules
        END DO

     ! End loop over all nitrogen molecules
  END DO


  ! Multiply energy and pressure by the L-J factor
  potenergy = 4.0d0* eps * potenergy
  virial = 8.0d0* eps * virial


END SUBROUTINE total_energy

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

SUBROUTINE adjust_mc_parameters

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Adjusts MC parameter maxdisp                                       !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Adjust maxdisp
  IF ( (DBLE(acctrans)/DBLE(ntrans)) > 0.50d0) THEN
     IF (maxdisp < box)  maxdisp = 1.050d0 * maxdisp
  ELSE
     maxdisp = 0.950d0 * maxdisp
  ENDIF

  RETURN

END SUBROUTINE adjust_mc_parameters

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  


SUBROUTINE sample_gr

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Updates the histogram for g(r)                                     !
  !                                                                    !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Local variables
  INTEGER :: imol,j
  DOUBLE PRECISION :: rxi,ryi,rzi
  DOUBLE PRECISION :: rxij,ryij,rzij,rijsq,rij
  INTEGER :: bin

  totalsamplesgr = totalsamplesgr + 1

  ! Start loop over all nitrogen molecules
  DO imol=1,n_molecules-1

     rxi = rx(imol)
     ryi = ry(imol)
     rzi = rz(imol)

     ! Contribution due to interactions with other nitrogen molecules

     ! Start loop over other molecules
     DO j=imol+1,n_molecules

        ! calculate vector between centers of mass
        rxij = rx(j) - rxi
        ryij = ry(j) - ryi
        rzij = rz(j) - rzi

        ! apply minimum image convention
        rxij = rxij - DNINT(rxij/box)*box
        ryij = ryij - DNINT(ryij/box)*box
        rzij = rzij - DNINT(rzij/box)*box

        ! Calculate the distance between centers of mass
        rijsq = rxij**2 + ryij**2 + rzij**2
        rij = dSQRT(rijsq)

        IF (rij < box/2) THEN
           bin = INT(rij/delr) + 1
           hist(bin) = hist(bin) + 2
        END IF

        ! End loop over other molecules
     END DO

     ! End loop over molecules
  END DO


  RETURN

END SUBROUTINE sample_gr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

SUBROUTINE start

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Sets everything up to start running the MC code                    !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

  USE global_variables
  IMPLICIT NONE

  ! Local variables
  INTEGER :: imol,bin
  CHARACTER(LEN=5) :: atype

  ! Read the input file
  OPEN(UNIT=10,FILE='inputfile.txt',STATUS='OLD',ACCESS='SEQUENTIAL',&
       &ACTION='READ')
  READ(10,*)
  READ(10,*) temp
  READ(10,*) cutoff_nn
  READ(10,*)
  READ(10,*)
  READ(10,*) sig, eps
  cutoff_nn_square = cutoff_nn * cutoff_nn
  READ(10,*)
  READ(10,*)
  READ(10,*) nblocks
  READ(10,*) nsteps
  READ(10,*) nadjust
  READ(10,*) nsamplegr
  READ(10,*) dcd_freq
  READ(10,*) maxdisp
  
  CLOSE(UNIT=10)


  ! Calculate the very small distances that will cause huge energies
  huge_energy = temp * huge_number
  small_rsq_nn = (0.80d0)**2

  ! Get initial positions of nitrogen molecules
  OPEN(UNIT=20,FILE='initial.xyz',STATUS='OLD',ACCESS='SEQUENTIAL',&
       &ACTION='READ')
  READ(20,*) n_molecules
  READ(20,*) box

  ! Read positions
  DO imol=1,n_molecules
     READ(20,*) atype,rx(imol),ry(imol),rz(imol)
  END DO
  CLOSE(UNIT=20)


  ! Initialize dummy variable for random number generator
  CALL initialize_idum(idum)

  ! Initialize acc/rej statistics
  acctrans = 0
  ntrans = 0

  ! Initialize radial distribution function
  DO bin=1,maxbin
     hist(bin) = 0
  END DO
  delr = box/(2*maxbin)
  constant=4*pi*DBLE(n_molecules)/(3.0d0*box**3)
  totalsamplesgr = 0

  ! Open the DCD file for writing
  CALL open_dcd

  ! Open file to write results
  OPEN(UNIT=40,FILE='blkaverages',STATUS='UNKNOWN',ACCESS='SEQUENTIAL')
  WRITE(40,*) '#NVT Results (blocks averages)'
  WRITE(40,*)
  WRITE(40,FMT='(F14.6, '' = Temperature'')') temp
  WRITE(40,FMT='(9X,I5, '' = Number of molecules'')') n_molecules
  WRITE(40,FMT='(F14.6, '' = Volume of the simulation box'')') box**3
  WRITE(40,*)
  IF(sig == 1.0 .AND. eps == 1.0) THEN
    WRITE(40,*) '#block       U (reduced)       U^2 (reduced)          P (reduced) '
  ELSE
    WRITE(40,*) 'block       U (J/mol)         U^2 (J/mol)^2          P (bar)'
  ENDIF
  
  RETURN

END SUBROUTINE start

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
赞一下 回复此楼
Life, Love, Laugh.
5楼2010-11-25 01:00:20
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
(再续前)



SUBROUTINE finish

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! Reports the final results of the MC code                           !
  !                                                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                    !
  ! qphll                                                              !
  !                                                                                                                                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  USE global_variables
  IMPLICIT NONE

  ! Local variables
  INTEGER :: imol,j
  DOUBLE PRECISION :: rlow,rup,nideal
  INTEGER :: bin


  OPEN(UNIT=50,FILE='final.xyz',STATUS='UNKNOWN',ACCESS='SEQUENTIAL')
  WRITE(50,FMT='(1I4,12X)') n_molecules
  WRITE(50,*) box,box,box
  DO imol=1,n_molecules
     WRITE(50,FMT='(A5,3F12.6)')'Ar', rx(imol),ry(imol),rz(imol)
  END DO

  CLOSE(UNIT=50)



  ! Calculate radial distribution function

  OPEN(UNIT=60,FILE='gr',STATUS='UNKNOWN',ACCESS='SEQUENTIAL')

  DO bin=1,maxbin
     rlow = delr*(bin-1)
     rup = rlow + delr
     nideal = constant*(rup**3-rlow**3)
     WRITE(60,FMT='(F8.4,F12.6)') rlow + delr/2 , &
          &hist(bin)/nideal/DBLE(n_molecules)/DBLE(totalsamplesgr)
  END DO

  CLOSE(UNIT=60)
     


  RETURN

END SUBROUTINE finish


  SUBROUTINE open_dcd

    USE global_variables
    IMPLICIT NONE


    ! Internal
    INTEGER :: j,n_frames
    INTEGER :: ntitle
    CHARACTER(LEN=80),DIMENSION(100) :: title
    INTEGER,DIMENSION(20) :: ictrl
    REAL*4 :: delta
    DOUBLE PRECISION,DIMENSION(6) :: pbcs

    ! Initialize delta
    delta = 1.0000

    ! Initialize ictrl
    ictrl = 0

    ! Calculate the number of frames
    n_frames = nblocks/dcd_freq

    ! Set ictrl elements
    ictrl(1) = n_frames
    ictrl(3) = 1
    ictrl(4) = ictrl(1)
    ictrl(11) = 0
    ictrl(20) = 24

    ! Set info for filler junk  
    ntitle = 2
    title(1) = "Created by JDM plugin                     "
    title(2) = "REMARKS Created 10 November, 2008 at 19:24"

    ! Set PBC information
    pbcs(1) = DBLE(box)
    pbcs(2) = 0.0d0
    pbcs(3) = DBLE(box)
    pbcs(4) = 0.0d0
    pbcs(5) = 0.0d0
    pbcs(6) = DBLE(box)

    OPEN(UNIT=5010,FILE='movie.dcd',FORM='UNFORMATTED',STATUS='UNKNOWN')

      ! Write the junk at the begining
      WRITE(5010) "CORD",(ictrl(j),j=1,9),delta, (ictrl(j),j=11,20)
      WRITE(5010) ntitle,(title(j),j=1,ntitle)
      WRITE(5010) n_molecules
  !   WRITE(5010) (pbcs(j),j=1,6)

    RETURN

  END SUBROUTINE open_dcd


  SUBROUTINE write_dcd(iblock)

    USE global_variables
    IMPLICIT NONE

    ! Passed
    INTEGER :: iblock

    ! Internal
    INTEGER :: j, n_frames, iframe

    ! Calculate the number of frames
    n_frames = nblocks/dcd_freq

    ! Calculate the current frame number
    iframe = iblock/dcd_freq

    ! Check if a step should be written at iblock
    IF(MOD(iblock,dcd_freq) == 0 ) THEN

      ! Write the coordinates to file (must be REAL*4)
      WRITE(5010)(REAL(rx(j)),j=1,n_molecules)
      WRITE(5010)(REAL(ry(j)),j=1,n_molecules)
      WRITE(5010)(REAL(rz(j)),j=1,n_molecules)

    END IF

    ! Close the file if it is the last frame to be written
    IF(iframe == n_frames) CLOSE(UNIT=5010)

    RETURN

  END SUBROUTINE write_dcd
赞一下 回复此楼
Life, Love, Laugh.
6楼2010-11-25 01:00:56
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
将前面三个中的内容合并存档为 qphll.f90 (哈哈, 必须保存为这个名字, 否则无法编译 ^-^)

然后随便找个F90的编译器, 比如


# ifort qphll.f90 -o qphll

或者 如果是PGI的话:

# pgf90 -o qphll qphll.for

应该就能得到可执行程序的.

不保证没有啥错误, 也不打算将我这边的.f90文件上传. 有些努力, 你总是要自己付出的.
赞一下 回复此楼
Life, Love, Laugh.
7楼2010-11-25 01:06:06
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
有枪了, 得有靶.

here you are.

initial.xyz 文件

32
4.3088693801 4.3088693801 4.3088693801
Ar   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
Ar   0.774800000000000E+00  0.774800000000000E+00  0.000000000000000E+00
Ar   0.000000000000000E+00  0.774800000000000E+00  0.774800000000000E+00
Ar   0.774800000000000E+00  0.000000000000000E+00  0.774800000000000E+00
Ar   0.154960000000000E+01  0.000000000000000E+00  0.000000000000000E+00
Ar   0.232440000000000E+01  0.774800000000000E+00  0.000000000000000E+00
Ar   0.154960000000000E+01  0.774800000000000E+00  0.774800000000000E+00
Ar   0.232440000000000E+01  0.000000000000000E+00  0.774800000000000E+00
Ar   0.000000000000000E+00  0.154960000000000E+01  0.000000000000000E+00
Ar   0.774800000000000E+00  0.232440000000000E+01  0.000000000000000E+00
Ar   0.000000000000000E+00  0.232440000000000E+01  0.774800000000000E+00
Ar   0.774800000000000E+00  0.154960000000000E+01  0.774800000000000E+00
Ar   0.154960000000000E+01  0.154960000000000E+01  0.000000000000000E+00
Ar   0.232440000000000E+01  0.232440000000000E+01  0.000000000000000E+00
Ar   0.154960000000000E+01  0.232440000000000E+01  0.774800000000000E+00
Ar   0.232440000000000E+01  0.154960000000000E+01  0.774800000000000E+00
Ar   0.000000000000000E+00  0.000000000000000E+00  0.154960000000000E+01
Ar   0.774800000000000E+00  0.774800000000000E+00  0.154960000000000E+01
Ar   0.000000000000000E+00  0.774800000000000E+00  0.232440000000000E+01
Ar   0.774800000000000E+00  0.000000000000000E+00  0.232440000000000E+01
Ar   0.154960000000000E+01  0.000000000000000E+00  0.154960000000000E+01
Ar   0.232440000000000E+01  0.774800000000000E+00  0.154960000000000E+01
Ar   0.154960000000000E+01  0.774800000000000E+00  0.232440000000000E+01
Ar   0.232440000000000E+01  0.000000000000000E+00  0.232440000000000E+01
Ar   0.000000000000000E+00  0.154960000000000E+01  0.154960000000000E+01
Ar   0.774800000000000E+00  0.232440000000000E+01  0.154960000000000E+01
Ar   0.000000000000000E+00  0.232440000000000E+01  0.232440000000000E+01
Ar   0.774800000000000E+00  0.154960000000000E+01  0.232440000000000E+01
Ar   0.154960000000000E+01  0.154960000000000E+01  0.154960000000000E+01
Ar   0.232440000000000E+01  0.232440000000000E+01  0.154960000000000E+01
Ar   0.154960000000000E+01  0.232440000000000E+01  0.232440000000000E+01
Ar   0.232440000000000E+01  0.154960000000000E+01  0.232440000000000E+01

(自己处理一下格式).

形如:

32
Box_x, Box_y, Box_z
Ar, x, y,z
Ar, x,y,z
.
.
.
.
.
.
赞一下 回复此楼
Life, Love, Laugh.
8楼2010-11-25 01:08:58
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
有枪有靶, 得有瞄准器.

here you are

inputfile.txt:



         PHYSICAL PARAMETERS
2.0                          Temperature (reduced)
5.0                          cutoff distance (reduced) (Must be <= (1/2)*box)

         LENNARD JONES PARAMETERS
1.0 1.0                      Sigma (reduced) Eps/k_b (reduced)

         MONTE CARLO RUN PARAMETERS
1000                         number of blocks
5000                         number of moves per block
1                            Adjust max displacement frequency (blocks)
1                            g(r) sample frequency (blocks)
1                            Write to dcd file frequency (blocks)
0.5                          Initial maximum displacement (box length)
赞一下 回复此楼
Life, Love, Laugh.
9楼2010-11-25 01:10:24
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qphll

金虫 (正式写手)
打靶回来告诉我成绩~~

洗洗睡吧.

明天继续.
赞一下 回复此楼
Life, Love, Laugh.
10楼2010-11-25 01:11:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 qphll 的主题更新
811/9返回列表123456下一页末页
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[找工作] 化学博士 +9 实名心动国民萌� 2024-12-17 16/800 2024-12-24 09:38 by maoder
[有机交流] 3-溴呋喃到甲氧基呋喃 100+4 dyl7838 2024-12-20 7/350 2024-12-24 08:49 by 091602
[职场人生] 国企、民企行业不同，如何决择？ +11 haihe623 2024-12-19 13/650 2024-12-24 08:49 by yiran909
[论文投稿] 挑数据是造假吗 +11 电话的建设 2024-12-21 11/550 2024-12-24 08:33 by 我乘着风
[考博] 2025申博 +3 梁艺艺 2024-12-21 4/200 2024-12-24 06:20 by firen2020
[文学芳草园] 生活的意义 +7 宋昙 2024-12-20 7/350 2024-12-24 05:22 by tfang
[硕博家园] 纠结退学ing +20 幽皮皮皮 2024-12-22 35/1750 2024-12-23 23:33 by 你看笑靥如花
[教师之家] 咨询一下，是不是教授之间比较忌讳谈论文造假，因为可能很多人都是造假 +9 akslis2024 2024-12-21 10/500 2024-12-23 20:46 by zeolitess
[硕博家园] 电池检测Sci 期刊推荐 +3 设置昵称不 2024-12-23 3/150 2024-12-23 20:40 by q74prt@q
[教师之家] 有两个省级人才帽子跳槽找工作，难度怎么样 +6 akslis2024 2024-12-23 6/300 2024-12-23 20:34 by 甄士隐207
[论文投稿] 三个手稿 10+5 s爱吃芒果 2024-12-22 9/450 2024-12-23 16:32 by 世界之源
[找工作] 水电解制氢方向求职 +6 kyle1999 2024-12-19 13/650 2024-12-23 10:02 by kyle1999
[硕博家园] 读研之后发现: 圈子变得越来越小了... +13 丁柯翔鸭 2024-12-19 18/900 2024-12-22 23:06 by 丁柯翔鸭
[功能材料] 请教大神，带隙是DFT计算出来的还是实验测出来的？ 10+3 大力2010 2024-12-21 6/300 2024-12-22 18:08 by Ayase_lumi
[考博] 华南理工大学 “新能源交叉创新团队--主动安全”课题组招收海外联合培养博士生 +4 hubble 2024-12-20 5/250 2024-12-21 16:16 by 那片叶落
[教师之家] 今年的科研绩效只奖励Q1了 +12 holypower 2024-12-20 17/850 2024-12-21 16:13 by willbuilder
[论文投稿] 论文校稿 50+5 whale_full 2024-12-17 10/500 2024-12-20 16:39 by 北京莱茵润色
[论文投稿] Transportation Research Part C 投稿模板 5+3 Yetzirah 2024-12-20 4/200 2024-12-20 16:32 by 北京莱茵润色
[论文投稿] 毕业论文疑虑 20+6 风吹荷叶煞 2024-12-17 17/850 2024-12-18 17:35 by holypower
[有机交流] 求这个反应的反应条件 +6 戒骄戒躁+ 2024-12-17 6/300 2024-12-18 11:43 by diaodeng
信息提示
关闭
请填处理意见
关闭