24小时热门版块排行榜

返回列表

当前主题已经存档。

ycsylvia

金虫 (初入文坛)

应助: 0 (幼儿园)
金币: 1654.5
帖子: 47
在线: 27小时
虫号: 130192
注册: 2005-12-11
专业: 凝聚态物性 II ：电子结构

[交流] 【求助】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)"

Eb
close(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, renergy
write(70, "(2x,e12.4, 'the imaginary part energy is ', E12.4)"

ke, kenergy
   close(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 ]

回复此楼