Ewald potential

      program ewald
      implicit none

      real*8 eta
      parameter(eta=0.8d0)     !тот самый параметр,от которого почему-то зависит потенциал

      real*8 eps
      parameter(eps=1.0d-9)

      integer N_max
      parameter(N_max=100)

      integer H_max
      parameter(H_max=30)

      integer k_max
      parameter(k_max=3*H_max**2)

      integer n,j
      integer H,k,m,n_e
      integer k1,k2,k3
      integer m1,m2,m3
      integer i1,i2,i3

      real*8 R_k(-H_max:H_max,-H_max:H_max,-H_max:H_max,3)
      real*8 e(N_max,4)    !   coordinates and values of original charges
      real*8 a1(3),a2(3),a3(3)

      real*8 g_m(-H_max:H_max,-H_max:H_max,-H_max:H_max,3)
      real*8 g_m2(-H_max:H_max,-H_max:H_max,-H_max:H_max)
      real*8 a2_a3(3),a3_a1(3),a1_a2(3)
      real*8 b1(3),b2(3),b3(3)
      real*8 Q_r
      real*8 pi

      real*8 w1,w2

      real*8 V_1(0:k_max,H_max)
      real*8 V_2(0:k_max,H_max)
      real*8 W_1(0:H_max)
      real*8 W_2(0:H_max)
      real*8 V

!==========================================================

      open(10,file='original_charges.inp',status='old')
      read(10,*)n_e
      do n=1,n_e
        read(10,*)j,(e(n,k),k=1,4)
      enddo
      close(10)

      pi=4*datan(1.0d0)

      a1(1)=0.0d0
      a1(2)=2.0d0
      a1(3)=2.0d0
      a2(1)=2.0d0
      a2(2)=0.0d0
      a2(3)=2.0d0
      a3(1)=2.0d0
      a3(2)=2.0d0
      a3(3)=0.0d0

!==========================================================

!   find direct lattice vectors R_k
      do k1=-H_max,H_max
        do k2=-H_max,H_max
          do k3=-H_max,H_max
            do j=1,3
              R_k(k1,k2,k3,j)=k1*a1(j)+k2*a2(j)+k3*a3(j)
            enddo
          enddo
        enddo
      enddo

!=================================================================================
!   find W_1(H)
!=================================================================================

      do H=1,H_max
        do k=0,k_max
          V_1(k,H)=0.0d0
        enddo
      enddo

      do H=1,H_max
        do k1=-H,H
          do k2=-H,H
            do k3=-H,H
              do j=1,n_e

                k=k1**2+k2**2+k3**2

                if (k.eq.0.0d0) then
                  V_1(k,H)=e(2,4)*erfc(eta*sqrt(e(2,1)**2+
     &                     e(2,2)**2+e(2,3)**2))
                else

      V_1(k,H)=V_1(k,H)+e(j,4)*erfc(eta*sqrt
     &         ((R_k(k1,k2,k3,1)+e(j,1))**2+
     &         (R_k(k1,k2,k3,2)+e(j,2))**2+
     &         (R_k(k1,k2,k3,3)+e(j,3))**2))/
     &         sqrt((R_k(k1,k2,k3,1)+e(j,1))**2+
     &         (R_k(k1,k2,k3,2)+e(j,2))**2+
     &         (R_k(k1,k2,k3,3)+e(j,3))**2)

                endif
              enddo
            enddo
          enddo
        enddo
      enddo

      do H=0,H_max
        W_1(H)=0.0d0
      enddo

      do H=1,H_max
        do k=0,k_max
          W_1(H)=W_1(H)+V_1(k,H)
        enddo
        if (abs(W_1(H)-W_1(H-1)).lt.eps) then
          w1=W_1(H-1)
          write(*,*)H-1,w1
          go to 1
        endif
      enddo

!===========================================================================

!   vector product [a2,a3]
1     do i1=1,3
        i2=i1+1
        if(i2.gt.3)i2=i2-3
        i3=i2+1
        if(i3.gt.3)i3=i3-3
        a2_a3(i1)=a2(i2)*a3(i3)-a2(i3)*a3(i2)
      enddo

!   vector product [a3,a1]
      do i1=1,3
        i2=i1+1
        if(i2.gt.3)i2=i2-3
        i3=i2+1
        if(i3.gt.3)i3=i3-3
        a3_a1(i1)=a3(i2)*a1(i3)-a3(i3)*a1(i2)
      enddo

!   vector product [a1,a2]
      do i1=1,3
        i2=i1+1
        if(i2.gt.3)i2=i2-3
        i3=i2+1
        if(i3.gt.3)i3=i3-3
        a1_a2(i1)=a1(i2)*a2(i3)-a1(i3)*a2(i2)
      enddo

!   find Q_r
      Q_r=0.0d0
      do j=1,3
        Q_r=Q_r+a1(j)*a2_a3(j)
      enddo

!   find b1,b2,b3
      do j=1,3
        b1(j)=a2_a3(j)/Q_r
      enddo

      do j=1,3
        b2(j)=a3_a1(j)/Q_r
      enddo

      do j=1,3
        b3(j)=a1_a2(j)/Q_r
      enddo

!   find reciptorical lattice vectors g_m
      do m1=-H_max,H_max
        do m2=-H_max,H_max
          do m3=-H_max,H_max
              do j=1,3
                g_m(m1,m2,m3,j)=2*pi*(m1*b1(j)+m2*b2(j)+m3*b3(j))
              enddo
          enddo
        enddo
      enddo

!   find g_m**2
      do m1=-H_max,H_max
        do m2=-H_max,H_max
          do m3=-H_max,H_max
            g_m2(m1,m2,m3)=0.0d0
          enddo
        enddo
      enddo

      do m1=-H_max,H_max
        do m2=-H_max,H_max
          do m3=-H_max,H_max
            do j=1,3
              g_m2(m1,m2,m3)=g_m2(m1,m2,m3)+
     &        g_m(m1,m2,m3,j)*g_m(m1,m2,m3,j)
            enddo
          enddo
        enddo
      enddo

!=================================================================================
!   find W_2(H)
!=================================================================================

      do H=1,H_max
        do m=0,k_max
          V_2(m,H)=0.0d0
        enddo
      enddo

      do H=1,H_max
        do m1=-H,H
          do m2=-H,H
            do m3=-H,H
              do j=1,n_e

              m=m1**2+m2**2+m3**2

              if (g_m2(m1,m2,m3).ne.0.0d0) then

              V_2(m,H)=V_2(m,H)+(4.0d0*pi/Q_r)*e(j,4)/g_m2(m1,m2,m3)*
     &        dcos((g_m(m1,m2,m3,1)*e(j,1))+(g_m(m1,m2,m3,2)*
     &        e(j,2))+(g_m(m1,m2,m3,3)*e(j,3)))*
     &        exp(-g_m2(m1,m2,m3)/(4.0d0*(eta)**2))

              endif

              enddo
            enddo
          enddo
        enddo
      enddo

      do H=0,H_max
        W_2(H)=0.0d0
      enddo

      do H=1,H_max
        do m=0,k_max
          W_2(H)=W_2(H)+V_2(m,H)
        enddo
        if (abs(W_2(H)-W_2(H-1)).lt.eps) then
          w2=W_2(H-1)
          write(*,*)H-1,w2
          go to 2
        endif
      enddo

!==============================================================================

2     V=w1+w2-2*eta*e(1,4)/sqrt(pi)
      write(*,*)V

      end