! Ewaldsum fortran 90 code ! 4/23/01 NAWH ! Input : T1x T1y T1z (Cartesian components of lattice vectors) ! : T2x T2y T2z ! : T3x T3y T3z ! : nions (Number of ions in unit cell) ! : q(i), tau(1:3,i), i=1,nion (charge and fractional position ! of ions) ! Note: ion positions = tau(1,i)*T1 + tau(2,i)*T2 + tau(3,i)*T3 ! : eps (error tolerance) IMPLICIT NONE real :: pi,eps, t1(3), t2(3), t3(3), g1(3),g2(3),g3(3),volcry, arg,x,gexp real :: eta, totalcharge, g1m, g2m, g3m, t1m, t2m, t3m,gcut, tmax ,ebsl,seta real :: tpi,glast2, con, con2 , cccc , ewald , v(3), w(3), rmag2 , prod integer :: nion, i, j, k, ng, nt , mmm1, mmm2, mmm3 , a , b real, allocatable :: q(:), tau(:,:) Write(6,*) ' Enter T1x, T1y, T1z in bohr units' Read(5,*) t1(1),t1(2),t1(3) Write(6,*) ' Enter T2x, T2y, T2z in bohr units' Read(5,*) t2(1),t2(2),t2(3) Write(6,*) ' Enter T3x, T3y, T3z in bohr units' Read(5,*) t3(1),t3(2),t3(3) pi = acos(-1.0) volcry = t1(1)*(t2(2)*t3(3)-t2(3)*t3(2)) + & t1(2)*(t2(3)*t3(1)-t2(1)*t3(3)) + & t1(3)*(t2(1)*t3(2)-t2(2)*t3(1)) g1(1) = 2 * pi * (t2(2)*t3(3)-t2(3)*t3(2))/volcry g1(2) = 2 * pi * (t2(3)*t3(1)-t2(1)*t3(3))/volcry g1(3) = 2 * pi * (t2(1)*t3(2)-t2(2)*t3(1))/volcry g2(1) = 2 * pi * (t3(2)*t1(3)-t3(3)*t1(2))/volcry g2(2) = 2 * pi * (t3(3)*t1(1)-t3(1)*t1(3))/volcry g2(3) = 2 * pi * (t3(1)*t1(2)-t3(2)*t1(1))/volcry g3(1) = 2 * pi * (t1(2)*t2(3)-t1(3)*t2(2))/volcry g3(2) = 2 * pi * (t1(3)*t2(1)-t1(1)*t2(3))/volcry g3(3) = 2 * pi * (t1(1)*t2(2)-t1(2)*t2(1))/volcry volcry = abs(volcry) t1m = SQRT(DOT_PRODUCT(t1,t1)) t2m = SQRT(DOT_PRODUCT(t2,t2)) t3m = SQRT(DOT_PRODUCT(t3,t3)) g1m = SQRT(DOT_PRODUCT(g1,g1)) g2m = SQRT(DOT_PRODUCT(g2,g2)) g3m = SQRT(DOT_PRODUCT(g3,g3)) Write(6,*) 'Input total number of ions in unit cell' Read(5,*) nion Allocate(q(nion),tau(3,nion)) Write(6,*) 'For each ion, input q and tau (in fractional coordinates of T)' do i=1,nion read(5,*) q(i),tau(1,i),tau(2,i),tau(3,i) enddo Write(6,*) ' Input Gcut for reciprocal lattice sum and error tolerance' Read(5,*) gcut , ebsl tpi=2*pi con=volcry/(4*pi) con2=(4*pi)/volcry glast2=gcut*gcut gexp=-alog(ebsl) eta=glast2/gexp Write(6,*) ' eta value for this calculation' , eta cccc=sqrt(eta/pi) x=0 do i=1,nion x=x+q(i)**2 enddo totalcharge = sum(q) Write(6,*) ' Total charge = ', totalcharge ewald=-cccc*x-4*pi*(totalcharge**2)/(volcry*eta) tmax=sqrt(2*gexp/eta) seta=sqrt(eta)/2 mmm1=tmax/t1m+1.5 mmm2=tmax/t2m+1.5 mmm3=tmax/t3m+1.5 Write (6,*) ' lattice summation indices -- ', mmm1,mmm2,mmm3 do a = 1,nion do b = 1,nion v(:) = (tau(1,a)-tau(1,b))*t1(:) + (tau(2,a)-tau(2,b))*t2(:) & + (tau(3,a)-tau(3,b))*t3(:) prod=q(a)*q(b) do i = -mmm1, mmm1 do j = -mmm2, mmm2 do k = -mmm3, mmm3 if ((a.ne.b).or.((abs(i)+abs(j)+abs(k)).ne.0)) then w(:) = v(:) + i*t1 + j*t2 + k*t3 rmag2 = sqrt(DOT_PRODUCT(w,w)) arg=rmag2*seta ewald = ewald + prod*erfc(arg)/rmag2 endif enddo enddo enddo enddo enddo mmm1=gcut/g1m+1.5 mmm2=gcut/g2m+1.5 mmm3=gcut/g3m+1.5 Write(6,*) ' Reciprocal lattice summation indices --', mmm1,mmm2,mmm3 do i = -mmm1, mmm1 do j = -mmm2, mmm2 do k = -mmm3, mmm3 if ((abs(i)+abs(j)+abs(k)).ne.0) then w(:) = i*g1(:) + j*g2(:) + k*g3(:) rmag2=DOT_PRODUCT(w,w) x=con2*exp(-rmag2/eta)/rmag2 do a = 1,nion do b = 1,nion v(:) = tau(:,a)-tau(:,b) prod = q(a)*q(b) arg=tpi*(i*v(1)+j*v(2)+k*v(3)) ewald=ewald + x*prod*cos(arg) enddo enddo endif enddo enddo enddo Write(6,*) ' Ewald energy in Rydbergs', ewald stop end