| 查看: 1582 | 回复: 7 | ||||
| 当前主题已经存档。 | ||||
ycsylvia金虫 (初入文坛)
|
[交流]
【求助】Ewald方法计算静电能 程序
|
|||
|
版主， 请不要删除哦，我在这里求助， 是因为分子动力学的好像常用Ewald方法计算静电能 若用该方法计算离子晶体的静电能，实空间只计算一个primitive cell， 傅立叶空间计算6×6×6的supercell，帮我看看程序这样写可以吗？ MODULE PARA IMPLICIT NONE INTEGER, PARAMETER :: KMAX=6 integer, parameter :: maxk=1000 REAL :: EFF=1.3899E-7 ! EFF=e*e*avgadelo constant/(4*pi*epsilon) Kj*m/mol REAL :: A=5.64E-10 !SI unit REAL :: EPSILON=1.0 REAL :: V=179.406E-30 !L*(A**3) integer, parameter :: ncell=6 integer, parameter :: nmolecular=7 integer, parameter :: nion=1512 !=nsublattice*nmolecule integer, parameter :: nsublattice=216 !=ncell**3 real :: qi, qj real ::rx(1:nion), ry(1:nion), rz(1:nion), q(1:nion) REAL :: KVEC(1000) END MODULE PROGRAM LRO_ALLOY USE PARA IMPLICIT NONE integer i,j,k,n real Eb real,external :: Hq,HqT n=1 do i=1, ncell do j=1, ncell do k=1, ncell !coordinates and charge of ions in sublattice rx(n)=real(i) ry(n)=real(j) rz(n)=real(k) q(n)=-6.0 n=n+1 rx(n)=real(i) ry(n)=real(j) rz(n)=real(k)+0.5 q(n)=+1.0 n=n+1 rx(n)=real(i) ry(n)=real(j) rz(n)=real(k)-0.5 q(n)=+1.0 n=n+1 rx(n)=real(i)+0.5 ry(n)=real(j) rz(n)=real(k) q(n)=+1.0 n=n+1 rx(n)=real(i)-0.5 ry(n)=real(j) rz(n)=real(k) q(n)=+1.0 n=n+1 rx(n)=real(i) ry(n)=real(j)+0.5 rz(n)=real(k) q(n)=+1.0 n=n+1 rx(n)=real(i) ry(n)=real(j)-0.5 rz(n)=real(k) q(n)=+1.0 n=n+1 end do end do end do n=n-1 Eb=0.0 Eb=Hq( ) open(60, file=" the electro energy.dat", status='replace') write(60, "('the total ion number is ', i5)" ) n write(60, "('the electrostatic energy is ', E12.4)"  Ebclose(60) STOP END FUNCTION Hq( ) USE PARA IMPLICIT NONE REAL KALPHA REAL RENERGY,KENERGY, Hq, re, ke INTEGER LX,LY,LZ, KN INTEGER MCSTEP REAL, EXTERNAL :: RWALD, KWALD integer time(0:1), timediff KALPHA = 1.0 RENERGY=0.0 KENERGY=0.0 Hq=0.0 re=RWALD(KALPHA) RENERGY=re*eff ke=KWALD(KALPHA) KENERGY=ke*eff Hq=(RENERGY+KENERGY)/nsublattice open(70, file=" the part energy.dat", status='replace') write(70, "(2x, e12.4,'the real part energy is ', E18.4)" re, renergywrite(70, "(2x,e12.4, 'the imaginary part energy is ', E12.4)" ke, kenergyclose(70) RETURN END SUBROUTINE SETUP ( KAPPA ) use para ! ******************************************************************* ! ** ROUTINE TO SET UP THE WAVE-VECTORS FOR THE EWALD SUM. ** ! ** ** ! ** THE WAVEVECTORS MUST FIT INTO A BOX OF UNIT LENGTH. ** ! ** IN THIS EXAMPLE WE ALLOW A MAXIMUM OF 1000 WAVEVECTORS. ** ! ******************************************************************* REAL KAPPA INTEGER KSQMAX, KSQ, KX, KY, KZ, TOTK REAL TWOPI, B, RKX, RKY, RKZ, RKSQ PARAMETER ( KSQMAX = 42 , TWOPI = 6.2831853 ) ! ******************************************************************* B = 1.0 / 4.0 / KAPPA / KAPPA ! ** LOOP OVER K-VECTORS. NOTE KX IS NON-NEGATIVE ** TOTK = 0 DO KX = 1, KMAX RKX = TWOPI * REAL ( KX ) DO KY = 1, KMAX RKY = TWOPI * REAL ( KY ) DO KZ = 1, KMAX RKZ = TWOPI * REAL ( KZ ) KSQ = KX * KX + KY * KY + KZ * KZ IF ( ( KSQ .LT. KSQMAX ) .AND. ( KSQ .NE. 0 ) ) THEN TOTK = TOTK + 1 IF ( TOTK .GT. maxk ) STOP 'KVEC IS TOO SMALL' RKSQ = RKX * RKX + RKY * RKY + RKZ * RKZ KVEC(TOTK) = TWOPI * EXP ( -B * RKSQ ) / RKSQ ENDIF ENDDO ENDDO ENDDO ! WRITE( *, ' ( '' EWALD SUM SETUP COMPLETE '' ) ' ) ! WRITE( *, ' ( '' NUMBER OF WAVEVECTORS IS '', I5 ) ' ) TOTK RETURN END subroutine FUNCTION RWALD ( KAPPA ) USE PARA ! ******************************************************************* ! ** CALCULATES R-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD. ** ! ** ** ! ** PRINCIPAL VARIABLES: ** ! ** ** ! ** INTEGER N NUMBER OF IONS ** ! ** REAL RX(N),RY(N),RZ(N) POSITIONS OF IONS ** ! ** REAL Z(N) IONIC CHARGES ** ! ** REAL VR R-SPACE POTENTIAL ENERGY ** ! ** REAL FUNCTION ERFC ( X ) ** ! ** RETURNS THE COMPLEMENTARY ERROR FUNCTION ** ! ******************************************************************* REAL KAPPA, VR REAL RXI, RYI, RZI, ZI, RXIJ, RYIJ, RZIJ REAL RIJSQ, RIJ, KRIJ, ERFC, VIJ ! INTEGER I, J ! ******************************************************************* DO I = 1, nion - 1 RXI = RX(I) RYI = RY(I) RZI = RZ(I) ZI = q(I) DO J = I + 1, nion RXIJ = RXI - RX(J) RYIJ = RYI - RY(J) RZIJ = RZI - RZ(J) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ RIJ = SQRT ( RIJSQ ) if (RIJ/=0.0) then KRIJ = KAPPA * RIJ VIJ = ZI * q(J) * ERFC ( KRIJ ) / RIJ VR = VR + VIJ endif ENDDO ENDDO RWALD=VR/A RETURN END FUNCTION KWALD( KAPPA ) USE PARA ******************************************************************* ! ** CALCULATES K-SPACE PART OF POTENTIAL ENERGY BY EWALD METHOD. ** ** ! ** INTEGER N NUMBER OF IONS ** ! ** REAL RX(N),RY(N),RZ(N) POSITIONS OF IONS ** ! ** REAL Z(N) IONIC CHARGES ** ! ** REAL VK K-SPACE POTENTIAL ENERGY ** ! ** REAL VKS SELF PART OF K-SPACE SUM ** ! ******************************************************************* ! REAL KVEC(MAXK), RX(N), RY(N), RZ(N), Z(N) REAL KWALD, KAPPA INTEGER TOTK INTEGER KX, KY, KZ, KSQMAX, KSQ,KXYZ REAL TWOPI, FACTOR, VD, VS, RSQPI PARAMETER ( KSQMAX = 42 ) PARAMETER ( TWOPI = 6.2831853, RSQPI = 0.5641896 ) COMPLEX EIKX(1:nion, 1:KMAX) COMPLEX EIKY(1:nion, 1:KMAX) COMPLEX EIKZ(1:nion, 1:KMAX) COMPLEX EIKR(nion), SUM ! ** CONSTRUCT EXP(IK.R) FOR ALL IONS AND K-VECTORS ** ! ** CALCULATE KX, KY, KZ = 0 , -1 AND 1 EXPLICITLY ** DO I = 1, nion EIKX(I, 1) = CMPLX ( COS ( TWOPI * RX(I) ) ,SIN ( TWOPI * RX(I) ) ) EIKY(I, 1) = CMPLX ( COS ( TWOPI * RY(I) ) ,SIN ( TWOPI * RY(I) ) ) EIKZ(I, 1) = CMPLX ( COS ( TWOPI * RZ(I) ) ,SIN ( TWOPI * RZ(I) ) ) ENDDO ! ** CALCULATE REMAINING KX, KY AND KZ BY RECURRENCE ** DO KX = 2, KMAX DO I = 1, nion EIKX(I, KX) = EIKX(I, KX-1) * EIKX(I, 1) ENDDO ENDDO DO KY = 2, KMAX DO I = 1, nion EIKY(I, KY) = EIKY(I, KY-1) * EIKY(I, 1) ENDDO ENDDO DO KZ = 2, KMAX DO I = 1, nion EIKZ(I, KZ) = EIKZ(I, KZ-1) * EIKZ(I, 1) ENDDO ENDDO ! ** SUM OVER ALL VECTORS ** VD = 0.0 TOTK = 0 DO KX = 1, KMAX DO KY = 1, KMAX DO KZ = 1, KMAX KSQ = KX * KX + KY * KY + KZ * KZ IF ( ( KSQ .LT. KSQMAX ) .AND. ( KSQ .NE. 0 ) ) THEN TOTK = TOTK + 1 SUM = (0.0, 0.0) DO I = 1, nion EIKR(I) = EIKX(I, KX) * EIKY(I, KY) * EIKZ(I, KZ) SUM = SUM + q(I) * EIKR(I) ENDDO VD = VD + KVEC(TOTK) * CONJG ( SUM ) * SUM ENDIF ENDDO ENDDO ENDDO ! ** CALCULATES SELF PART OF K-SPACE SUM ** VS = 0.0 DO I = 1, nion VS = VS + q(I) * q(I) ENDDO VS = RSQPI * KAPPA * VS !vs=sqrt(alpha/pi)*qi*qi sum over all i ! ** CALCULATE THE TOTAL K-SPACE POTENTIAL ** KWALD = VD/v- VS RETURN END real function erfc(x) ! # MS Fortran ! Complementary error function from Numerical Recipes. implicit none real :: x real :: t,z z = abs(x) t = 1.0 / ( 1.0 + 0.5 * z ) erfc = t * exp( -z * z - 1.26551223 + t * & ( 1.00002368 + t * ( 0.37409196 + t * & ( 0.09678418 + t * (-0.18628806 + t * & ( 0.27886807 + t * (-1.13520398 + t * & ( 1.48851587 + t * (-0.82215223 + t * 0.17087277 ))))))))) if ( x.lt.0.0 ) erfc = 2.0 - erfc end function erfc [ Last edited by lei0736 on 2009-11-25 at 12:28 ] |
回复此楼» 收录本帖的淘帖专辑推荐
 分子动力学 |
» 猜你喜欢
[复试调剂]西南科技大学国防/材料导师推荐
已经有6人回复
-
化学工程321分求调剂
已经有12人回复
-
211本，11408一志愿中科院277分，曾在中科院自动化所实习
已经有4人回复
-
材料专硕326求调剂
已经有5人回复
-
东南大学364求调剂
已经有5人回复
-
国自科面上基金字体
已经有7人回复
-
药学383 求调剂
已经有4人回复
-
286求调剂
已经有5人回复
-
085601求调剂
已经有3人回复
-
302求调剂
已经有5人回复
» 本主题相关商家推荐: (我也要在这里推广)
2楼2009-05-12 08:41:16
superdirac
木虫 (正式写手)
- 应助: 2 (幼儿园)
- 金币: 4269.5
- 散金: 33
- 红花: 2
- 帖子: 536
- 在线: 287.5小时
- 虫号: 522969
- 注册: 2008-03-12
- 专业: 理论和计算化学
3楼2009-05-12 10:01:50
ycsylvia
金虫 (初入文坛)
- 应助: 0 (幼儿园)
- 金币: 1654.5
- 帖子: 47
- 在线: 27小时
- 虫号: 130192
- 注册: 2005-12-11
- 专业: 凝聚态物性 II ：电子结构
4楼2009-05-12 15:18:46
softmatter
木虫 (正式写手)
- 应助: 2 (幼儿园)
- 金币: 4439.2
- 帖子: 723
- 在线: 73.5小时
- 虫号: 298219
- 注册: 2006-11-19
- 专业: 凝聚态物性I:结构、力学和
5楼2009-05-12 17:16:09
6楼2009-05-12 18:23:37
7楼2009-07-21 11:37:04
bay__gulf
金虫 (著名写手)
刘苏州
- 模拟EPI: 8
- 应助: 9 (幼儿园)
- 贵宾: 4.85
- 金币: 2332.8
- 红花: 1
- 帖子: 1344
- 在线: 271小时
- 虫号: 592012
- 注册: 2008-09-03
- 专业: 理论和计算化学
- 管辖: 分子模拟
8楼2009-07-21 12:43:39
