| 查看: 8716 | 回复: 80 | ||||||||||||
| 本帖产生 2 个 模拟EPI ，点击这里进行查看 | ||||||||||||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||||||||||||
qphll金虫 (正式写手)
|
[交流]
【原创】抛土引玉, Monte Carlo 如何自己写程序: NVT_LJ已有57人参与
|
|||||||||||
|
引子: 注册小木虫很多年了, 一直是条小虫子, 也喜欢这样的论坛, 有虫, 有大虫, 有菜, 也有肉. 在学术和科研这个大盘子上面, 谁都是小虫子, 谁都能成为某一方向的专家.或许大家内心偶尔还是有那种"为中华之崛起而读书"的热血冲动, 但是更多时候, 为现实的残酷而黯然, 为学海之无涯而沮丧. 但是正像虫卵可以慢慢长成小虫一样, 只要你不断努力, 小虫也会破茧成为美丽的蝴蝶. 不管哭与笑, 不论败与成, Keep fighting, Never give up. 与大家共勉. qphll Nov. 24, 2010 |
回复此楼» 收录本帖的淘帖专辑推荐
 材料计算模拟实用技巧 | 分子模拟 | 生活休闲 | 学习方法 |
| condensed matter physics | 第一性原理和电化学 | 量子化学理论及量化软件专辑 | 量化,第一性,MD笔记 |
| 分子生物物理理论 |
» 本帖已获得的红花（最新10朵）
honglyn» 猜你喜欢
离职
已经有6人回复
-
求老师哦！
已经有9人回复
-
要不要加博导微信
已经有9人回复
-
四线城市公办高校无编or二线城市民办高校
已经有35人回复
-
你们都降薪了吗？
已经有16人回复
-
求助，要不要根据审稿人的要求引用他们的文献？
已经有7人回复
-
纠结退学ing
已经有37人回复
-
自己找博士可真难啊
已经有16人回复
-
只想摆烂
已经有19人回复
-
研一新生（无科研经历）求助
已经有5人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- Handbook of Monte Carlo Methods 蒙托卡罗法手册
已经有159人回复
- Monte Carlo模拟高分子玻璃化转变
已经有10人回复
- fortran call语句
已经有5人回复
- 关于fix nvt里的Tdamp
已经有7人回复
- 分子动力学和Monte Carlo模拟-南开大学讲义
已经有121人回复
- Monte carlo模拟
已经有5人回复
- 关于Reverse Monte Carlo---2
已经有7人回复
- RAN(IDUM)源代码
已经有3人回复
- Monte Carlo 模拟时图是用什么软件做的
已经有10人回复
- 国内有哪些单位做用Monte Carlo模拟催化反应过程的
已经有3人回复
- Numerical Recipes in Fortran
已经有181人回复
- 关于 “Just Accepted”
已经有5人回复
- 【求助】是不是monte carlo 模拟只适用于高分子方向的？
已经有6人回复
- 【求助】求如何学习monte carlo仿真
已经有7人回复
- 【求助】用monte carlo法做分子模拟，用什么编程语言实现比较好？
已经有25人回复
- 弱弱的问句：权威，一级，二级，CSSCI如何划分？
已经有9人回复
- 【求助】Total energy 是什么的能量
已经有5人回复
- 【转帖】2维无外场Ising模型的Monte Carlo模拟的C++实现
已经有7人回复
- 【求助】想用Monte Carlo方法预测晶体磁性转变温度应该用什么程序
已经有5人回复
- 【求助】cutoff distance
已经有4人回复
- 【在线答疑】经典粒子体系的Monte Carlo 模拟之基础篇
已经有73人回复
qphll
金虫 (正式写手)
- 模拟EPI: 10
- 应助: 18 (小学生)
- 金币: 2708.7
- 散金: 2294
- 红花: 8
- 帖子: 554
- 在线: 227.1小时
- 虫号: 89654
- 注册: 2005-08-29
- 性别: GG
- 专业: 碳素材料与超硬材料
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10, 模拟EPI+1):鼓励 2010-11-25 09:24:33
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 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 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 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 |
4楼2010-11-25 00:58:35
10