!  Copyright (C) 2017 The authors of Py-ChemShell
!
!  This file is part of Py-ChemShell.
!
!  Py-ChemShell is free software: you can redistribute it and/or modify
!  it under the terms of the GNU Lesser General Public License as
!  published by the Free Software Foundation, either version 3 of the
!  License, or (at your option) any later version.
!
!  Py-ChemShell is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with Py-ChemShell.  If not, see
!  <http://www.gnu.org/licenses/>.

subroutine periodic_electrostatics(grid_mode, ncores, nbqs, nshells, ngridpts, ierr)

  ! Calculate the electrostatic potential 
  ! and electric field for a periodic system's unit cell and store in field object
  
  ! grid_mode = 0: at the core and shell positions
  ! grid_mode = 1: at the positions specified by a grid field
  
  use ChemShellClusterModule, only : get_ndims, debug
  
  implicit none

  integer, intent(in) :: grid_mode, ncores, nbqs, nshells, ngridpts
  integer, intent(out) :: ierr

  integer :: ndim, nresults, soffset, send, i

  double precision, dimension(3, ncores+nshells) :: cart_coords
  double precision, dimension(3, ngridpts) :: cart_grid

  double precision, allocatable :: pot(:)
  double precision, allocatable :: efield(:,:)

  integer, parameter :: periodic_frag = 1, cluster_frag = 2

  integer, external :: set_grid_electrostatics

  ierr = 0

  ndim = get_ndims()

  ! Allocate results arrays and get grid coordinates if necessary
  if(debug.gt.0) write(6,*) 'Grid mode set to', grid_mode
  select case (grid_mode)
  case (0)
     cart_grid = 0.0d0
     nresults = ncores + nshells
  case (1)
     call get_grid_coords(ngridpts, cart_grid, ierr)
     if (ierr .ne. 0) return
     nresults = ngridpts
  end select
  allocate(pot(nresults))
  allocate(efield(3,nresults))

  if(debug.gt.0) write(6,*) 'Number of points for FIELD object=', nresults

  call periodic_electrostatics_f(grid_mode, ndim, ncores, nbqs, nshells, ngridpts, &
       cart_grid, nresults, pot, efield, ierr)
  if (ierr .ne. 0) return

  ! Write the potential and electric field to field objects
  select case (grid_mode)
  case (0)
     call get_cart_coords(periodic_frag, cart_coords, ncores, nshells, ierr)
     if (ierr .ne. 0) return
     soffset = ncores + 1
     send = ncores + nshells
     if(debug.gt.0) then
       print *
       print *, '--------------------------------'
       print *, 'Values of potentials:'
       do i = 1, nresults
         print *, pot(i)
       end do
       print *, '--------------------------------'
       print *
       print *, 'values of efield:',efield(:,:)
     end if
     call set_electrostatics(periodic_frag, cart_coords, pot(1:ncores), &
          pot(soffset:send), efield(1:3,1:ncores), efield(1:3,soffset:send), &
          ncores, nshells, ierr)
     if (ierr .ne. 0) return
  case (1)
     if(set_grid_electrostatics(ngridpts, cart_grid, pot, efield) .ne. 0) then
        write(6, *) ">>> ChemShell CLUSTER ERROR: set_grid_electrostatics called in periodic_electrostatics"
     end if
  end select

  deallocate(pot)
  deallocate(efield)

end subroutine periodic_electrostatics


subroutine periodic_electrostatics_f(grid_mode, ndim, ncores, nbqs, nshells, ngridpts, &
     cart_grid, nresults, pot, efield, ierr)
  !
  ! Calculate the electrostatic potential and electric field for a periodic system's
  ! unit cell
  !
  ! grid_mode = 0: at the core and shell positions
  ! grid_mode = 1: at the positions specified by supplied grid coordinates
  !
  ! Can be called directly from Fortran routines to avoid faffing around with field objects
  ! 
  implicit none

  integer, intent(in) :: grid_mode, ndim, ncores, nbqs, nshells, ngridpts, nresults
  double precision, dimension(3, ngridpts), intent(in) :: cart_grid
  double precision, dimension(nresults), intent(out) :: pot
  double precision, dimension(3,nresults), intent(out) :: efield
  integer, intent(out) :: ierr

  integer :: i, j, k, ishel, pshel, ashel, isite

  double precision :: a, b, c, alpha, beta, gamma, area, volume
  double precision, dimension(3, 3) :: vect, vecinv, cart_to_frac
  double precision, dimension(3, ncores+nshells) :: cart_coords
  double precision, dimension(3, ncores) :: f_core
  double precision, dimension(3, nshells) :: f_shel
  double precision, dimension(ncores+nshells) :: cart_charges
  double precision, dimension(ncores) :: charge_core
  double precision, dimension(nshells) :: charge_shel
  integer, dimension(nshells) :: point_core
  integer, dimension(ncores) :: point_shel
  double precision, dimension(3, ngridpts) :: f_grid

  double precision :: dum
  double precision, dimension(2) :: dum2
  double precision, dimension(3) :: dum3

  logical :: shell_model

  double precision :: two_pi

  integer, parameter :: periodic_frag = 1, cluster_frag = 2

  ierr = 0

  ! PGI does not like acos in parameter initialisation so initialise here
  two_pi = 2.d0*dacos(-1.d0)

  ! Check there are no point charges
  if (nbqs .ne. 0) then
     write(6,*) "periodic_electrostatics_f: number of point charges must be zero"
     ierr = 1
     return
  end if

  ! Check nresults is consistent with grid_mode
  if ((grid_mode == 0 .and. nresults /= ncores + nshells) .or. &
       (grid_mode == 1 .and. nresults /= ngridpts)) then
     write(6,*) "periodic_electrostatics_f: inconsistent nresults"
     ierr = 1
     return
  end if
     
  if (nshells .eq. 0) then
     shell_model = .false.
  else
     shell_model = .true.
  end if

  ! Get unit cell vectors from fragment
  vect = 0.0d0
  call get_cell_vectors(periodic_frag, vect, ierr)
  if (ierr .ne. 0) return

  ! Call construct routines to generate unit cell parameters
  ! as they will not necessarily be present in the fragment object
  select case (ndim)
  case (2)
     call lattice_parameters_2(vect(1:2, 1:2), area, a, b, gamma)
  case (3)
     call lattice_parameters_3(vect, volume, a, b, c, alpha, beta, gamma)
  case default
     write(6,*) "periodic_electrostatics: error, ndim=", ndim
     ierr = 1
     return
  end select

  ! Calculate inverse cell vector
  vecinv = 0.0d0
  select case (ndim)
  case (2)
     call inverse_2(vect(1:2, 1:2), vecinv(1:2, 1:2), dum)
  case (3)
     call inverse_3(vect, vecinv)
  case default
     write(6,*) "periodic_electrostatics: error, ndim=", ndim
     ierr = 1
     return
  end select

  ! The construct routine takes fractional coordinates of the unit cell
  ! as an input. ChemShell fractional coordinates refer only to the asymmetric
  ! unit. They also are not 'the closest to the origin' (see frac_to_cart).
  ! Therefore we take the final cartesian coordinates and convert back
  ! to fractional coordinates.

  ! Get cartesian coordinates from fragment
  call get_cart_coords(periodic_frag, cart_coords, ncores, nshells, ierr)  
  if (ierr .ne. 0) return

  ! Convert from cartesians to fractional coordinates
  select case (ndim)
  case (2)
     do i = 1, 2
        do j = 1, 2
           cart_to_frac(i,j) = vecinv(j,i)/two_pi
        end do
     end do
     do isite = 1, ncores
        f_core(1, isite) = cart_to_frac(1,1) * cart_coords(1, isite) &
             + cart_to_frac(2,1) * cart_coords(2, isite)
        f_core(2, isite) = cart_to_frac(1,2) * cart_coords(1, isite) &
             + cart_to_frac(2,2) * cart_coords(2, isite)
        ! Just store the z component
        f_core(3, isite) = cart_coords(3, isite)
     end do
     do isite = 1, nshells
        ishel = ncores + isite
        f_shel(1, isite) = cart_to_frac(1,1) * cart_coords(1, ishel) &
             + cart_to_frac(2,1) * cart_coords(2, ishel)
        f_shel(2, isite) = cart_to_frac(1,2) * cart_coords(1, ishel) &
             + cart_to_frac(2,2) * cart_coords(2, ishel)
        ! Just store the z component
        f_shel(3, isite) = cart_coords(3, ishel)
     end do
     if (grid_mode == 1) then
        do isite = 1, ngridpts
           f_grid(1, isite) = cart_to_frac(1,1) * cart_grid(1, isite) &
                + cart_to_frac(2,1) * cart_grid(2, isite)
           f_grid(2, isite) = cart_to_frac(1,2) * cart_grid(1, isite) &
                + cart_to_frac(2,2) * cart_grid(2, isite)
           ! Just store the z component
           f_grid(3, isite) = cart_grid(3, isite)
        end do
     end if
  case (3)
     do i = 1, 3
        do j = 1, 3
           cart_to_frac(i,j) = vecinv(j,i)/two_pi
        end do
     end do
     do isite = 1, ncores
        f_core(1, isite) = cart_to_frac(1,1) * cart_coords(1, isite) &
             + cart_to_frac(2,1) * cart_coords(2, isite) + cart_to_frac(3,1) * cart_coords(3, isite)
        f_core(2, isite) = cart_to_frac(1,2) * cart_coords(1, isite) &
             + cart_to_frac(2,2) * cart_coords(2, isite) + cart_to_frac(3,2) * cart_coords(3, isite)
        f_core(3, isite) = cart_to_frac(1,3) * cart_coords(1, isite) &
             + cart_to_frac(2,3) * cart_coords(2, isite) + cart_to_frac(3,3) * cart_coords(3, isite)
     end do
     do isite = 1, nshells
        ishel = ncores + isite
        f_shel(1, isite) = cart_to_frac(1,1) * cart_coords(1, ishel) &
             + cart_to_frac(2,1) * cart_coords(2, ishel) + cart_to_frac(3,1) * cart_coords(3, ishel)
        f_shel(2, isite) = cart_to_frac(1,2) * cart_coords(1, ishel) &
             + cart_to_frac(2,2) * cart_coords(2, ishel) + cart_to_frac(3,2) * cart_coords(3, ishel)
        f_shel(3, isite) = cart_to_frac(1,3) * cart_coords(1, ishel) &
             + cart_to_frac(2,3) * cart_coords(2, ishel) + cart_to_frac(3,3) * cart_coords(3, ishel)
     end do
     if (grid_mode == 1) then
        do isite = 1, ngridpts
           f_grid(1, isite) = cart_to_frac(1,1) * cart_grid(1, isite) &
                + cart_to_frac(2,1) * cart_grid(2, isite) + cart_to_frac(3,1) * cart_grid(3, isite)
           f_grid(2, isite) = cart_to_frac(1,2) * cart_grid(1, isite) &
                + cart_to_frac(2,2) * cart_grid(2, isite) + cart_to_frac(3,2) * cart_grid(3, isite)
           f_grid(3, isite) = cart_to_frac(1,3) * cart_grid(1, isite) &
                + cart_to_frac(2,3) * cart_grid(2, isite) + cart_to_frac(3,3) * cart_grid(3, isite)
        end do
     end if
  case default
     write(6,*) "periodic_electrostatics: error, ndim=", ndim
     ierr = 1
     return
  end select

  ! Get charges of cartesian sites from fragment
  ! Split into core and shell arrays
  call get_cart_charges(periodic_frag, cart_charges, ncores+nshells, ierr)
  if (ierr .ne. 0) return
  charge_core = cart_charges(1:ncores)
  if (shell_model) charge_shel = cart_charges(ncores+1:ncores+nshells)

  ! Get pointers from shells to cores
  ! Then setup pointers in the other direction
  point_shel = 0
  if (shell_model) then
     call get_point_core(periodic_frag, point_core, nshells, ierr)
     if (ierr .ne. 0) return
     do i = 1, nshells
        point_shel(point_core(i)) = i
     end do
  end if

  select case (grid_mode)
  case (0)
     ! Calculate potential and field for all atoms and shells
     select case (ndim)
     case (2)
        dum = 0.0d0
        dum2 = 0.0d0
        do i = 1, ncores
           call periodic_pot_2d(dum2, dum, i, .false., &
                shell_model, ncores, nshells, point_shel, &
                f_core, f_shel, charge_core, charge_shel, &
                a, b, gamma, &
                vect(1:2, 1:2), vecinv(1:2, 1:2), area, &
                pot(i), efield(:, i))
           efield(:, i) = -efield(:, i)
           if (shell_model) then
              pshel = point_shel(i)
              if (pshel .ne. 0) then
                 ashel = pshel + ncores
                 call periodic_pot_2d(dum2, dum, i, .true., &
                      shell_model, ncores, nshells, point_shel, &
                      f_core, f_shel, charge_core, charge_shel, &
                      a, b, gamma, &
                      vect(1:2, 1:2), vecinv(1:2, 1:2), area, &
                      pot(ashel), efield(:, ashel))
                 efield(:, ashel) = -efield(:, ashel)
              end if
           end if
        end do
     case (3)
        dum3 = 0.0d0
        do i = 1, ncores
           call periodic_pot_3d(dum3, i, .false., &
                shell_model, ncores, nshells, point_shel, &
                f_core, f_shel, charge_core, charge_shel, &
                a, b, c, alpha, beta, gamma, &
                vect, vecinv, volume, &
                pot(i), efield(:, i))
           efield(:, i) = -efield(:, i)
           if (shell_model) then
              pshel = point_shel(i)
              if (pshel .ne. 0) then 
                 ashel = pshel + ncores
                 call periodic_pot_3d(dum3, i, .true., &
                      shell_model, ncores, nshells, point_shel, &
                      f_core, f_shel, charge_core, charge_shel, &
                      a, b, c, alpha, beta, gamma, &  
                      vect, vecinv, volume, &
                      pot(ashel), efield(:, ashel)) 
                 efield(:, ashel) = -efield(:, ashel)
              end if
           end if
        end do
     case default
        write(6,*) "periodic_electrostatics: error, ndim=", ndim
        ierr = 1
        return
     end select
  case (1)
     ! Calculate potential and field at grid points
     select case (ndim)
     case (2)
        do i = 1, ngridpts
           call periodic_pot_2d(f_grid(1:2, i), f_grid(3, i), 0, .false., &
                shell_model, ncores, nshells, point_shel, &
                f_core, f_shel, charge_core, charge_shel, &
                a, b, gamma, &
                vect(1:2, 1:2), vecinv(1:2, 1:2), area, &
                pot(i), efield(:, i))
           efield(:, i) = -efield(:, i)
        end do
     case (3)
        do i = 1, ngridpts
           call periodic_pot_3d(f_grid(:, i), 0, .false., &
                shell_model, ncores, nshells, point_shel, &
                f_core, f_shel, charge_core, charge_shel, &
                a, b, c, alpha, beta, gamma, &
                vect, vecinv, volume, &
                pot(i), efield(:, i))
           efield(:, i) = -efield(:, i)
        end do
     case default
        write(6,*) "periodic_electrostatics: error, ndim=", ndim
        ierr = 1
        return
     end select
  end select

end subroutine periodic_electrostatics_f


subroutine periodic_pot_2d(f_probe,z_probe,j_site,calculate_shell,&
     shell_model, no_core, no_shel, point_shel, &
     f_core, f_shel, charge_core, charge_shel, &
     a_cell, b_cell, gamma_cell, &
     vector, vecinv, area, &
     potential, efield)

  use ChemShellClusterModule, only : get_ndims, debug
!
! Calculates electrostatic potential at point f_probe, z_probe 
! if j_site is not 0, potential and field are calculated
! at a lattice site number j_site
!
  implicit none

  external dconerfc

  integer :: accuracy
  parameter (accuracy=20)
  double precision :: float_nr_nd
  parameter (float_nr_nd=1.1d0)
  ! Cutoff for exp*erfc term
  double precision :: xbig
  parameter (xbig = 26.543d0)

  logical :: shell_model, calculate_shell
  double precision, dimension(2) :: f_probe
  double precision :: z_probe
  double precision :: potential, potential_real_0, potential_real_1, &
       potential_real_2, potential_reciprocal
  double precision, dimension(3) :: efield, efield_real_1, &
       efield_real_2, efield_reciprocal
! NB: cartesian z-coordinate stored in f_...(3,:)
! cf construct, where it is passed in via r_...(3,:) instead
  double precision, dimension(3, no_core) :: f_core
  double precision, dimension(3, no_shel) :: f_shel
  double precision, dimension(2, no_core) :: g_core
  double precision, dimension(2, no_shel) :: g_shel
  double precision, dimension(no_core) :: charge_core
  double precision, dimension(no_shel) :: charge_shel

  double precision :: a_cell, b_cell, gamma_cell
  double precision, dimension(2, 2) :: vector, vecinv
  double precision :: area

  double precision :: pi

  double precision :: ewald_const, ewald_const_2, accuracy_factor, &
       rec_const, two_const, r_ewald_cut, r_v

  double precision :: dconerfc, fconerfc, fconerfc_1

  double precision :: distance_2, dist, ew_dist
  double precision, dimension(2) :: f, r, f_q, r_q, f_r, r_r
  double precision :: z, z_q, r_q_abs, z_abs, z_r
  double precision :: r_q_2, qr_c, qr_s, func_1, func_2, func_3, func_4, ftmp
  double precision :: contr_c, contr_s, factor

  double precision :: arg1, arg2, arg3, ha_to_ev, force_convert

  integer :: j_site, no_core, no_shel
  integer, dimension(no_core) :: point_shel
  integer :: max_r_x, max_r_y
  integer :: max_q_x, max_q_y
  integer :: ch_stdout, ix, iy
  integer :: j, k, n

  data ch_stdout/6/

  ! Choice of extents of the supercell in reciprocal space.
  ! Seems to be independent of direct latice with aproproate choice of Ewald's
  ! constant. The current choice provides sufficient precision for real*8.
  data max_q_x/10/, max_q_y/10/

  ! Initialise
  pi = dacos(-1.d0)
  ha_to_ev=27.211654d0
  force_convert=51.42206048632d0
  potential=0.d0
  potential_real_1=0.d0
  potential_real_2=0.d0
  potential_reciprocal=0.d0
  do k=1,3
     efield(k)=0.d0
     efield_real_1(k)=0.d0
     efield_real_2(k)=0.d0
     efield_reciprocal(k)=0.d0
  enddo
  ! Choose Ewald constant and set other constants
  ewald_const_2=float_nr_nd*pi/area
  ewald_const=sqrt(ewald_const_2)
  !
  ! Factor for reciprocal term is increased by a factor of 2 
  ! on collapsing complex exponents to cosine
  !
  rec_const=2.d0*pi/area
  two_const=2.d0*ewald_const/sqrt(pi)
  !
  ! Calculate extensions of UC, which includes circle of radius 
  ! sqrt(-ln(10^accuracy))/ewald_const
  ! 
  accuracy_factor=10.d0**accuracy
  accuracy_factor=sqrt(log(accuracy_factor))
  r_ewald_cut=accuracy_factor/ewald_const
  r_v=r_ewald_cut/abs(sin(gamma_cell))
  max_r_x=r_v/a_cell+1
  max_r_y=r_v/b_cell+1

  if(debug > 0 .and. j_site.eq.1) then
     write(*,'(a,g16.8,a,g16.8,a,2i4)') &
          'Ewald constant: ',ewald_const_2*3.57106767, &
          '  Ewald cut-off: ',r_ewald_cut*0.529177d0, &
          ' max: ',max_r_x,max_r_y, &
          'q sum factor:         ',rec_const
  endif

  if(j_site.ne.0) then
     n=point_shel(j_site)
     if(calculate_shell) then
        do k=1,2
           f_probe(k)=f_shel(k,n)
        enddo
        z_probe=f_shel(3,n)
     else
        do k=1,2
           f_probe(k)=f_core(k,j_site)
        enddo
        z_probe=f_core(3,j_site)
     endif
  endif
  ! Centre cell around f_probe
  ! with fractional coordinates named g_core & g_shel
  do j=1,no_core
     do k=1,2
        f(k)=f_core(k,j)
     enddo
     z=f_core(3,j)
     call translate_2(vector,f_probe,f,r)
     do k=1,2
        g_core(k,j)=f(k)
     enddo
     if(shell_model) then
        n=point_shel(j)
        if(n.ne.0) then
           do k=1,2
              g_shel(k,n)=f_shel(k,n)+(g_core(k,j)-f_core(k,j))
           enddo
        endif
     endif
  enddo
  ! 1st uc contribution in real space
  do j=1,no_core
     if(j.ne.j_site) then
        do k=1,2
           f(k)=g_core(k,j)
           f_r(k)=f_probe(k)-f(k)
        enddo
        z=f_core(3,j)
        z_r=z_probe-z
        call transform_2(vector,f_r,r_r)
        dist=distance_2(vector,f_probe,z_probe,f,z)
        ! this check is relevant for very large UC
        if(dist.lt.r_ewald_cut) then
           ew_dist=ewald_const*dist

           if (debug > 1 .and. j_site == 1) then
              write(*,'(3x,i4,a,5g20.12)') &
                   j,' c ',charge_core(j),dist*0.529177d0, &
                   (f_r(k),k=1,2),z_r*0.529177d0
           endif
           
           func_1=dconerfc(ew_dist)/dist
           potential_real_1=potential_real_1+ &
                charge_core(j)*func_1
           func_2=(two_const*exp(-ew_dist*ew_dist)+func_1)* &
                charge_core(j)/(dist*dist)
           do k=1,2
              efield_real_1(k)=efield_real_1(k)-func_2*r_r(k)
           enddo
           efield_real_1(3)=efield_real_1(3)-func_2*z_r
           if(shell_model) then
              n=point_shel(j)
              if(n.ne.0) then
                 do k=1,2
                    f(k)=g_shel(k,n)
                    f_r(k)=f_probe(k)-f(k)
                 enddo
                 z=f_shel(3,n)
                 z_r=z_probe-z
                 call transform_2(vector,f_r,r_r)
                 dist=distance_2(vector,f_probe,z_probe,f,z)

                 if (debug > 1) then
                    write(*,'(3x,i4,a,5f16.10)') &
                         j,' s ',charge_shel(n),dist*0.529177d0, &
                         (f_r(k),k=1,2),z_r*0.529177d0
                 end if

                 ew_dist=ewald_const*dist
                 func_1=dconerfc(ew_dist)/dist
                 potential_real_1=potential_real_1+ &
                      charge_shel(n)*func_1
                 func_2=(two_const*exp(-ew_dist*ew_dist)+func_1)* &
                      charge_shel(n)/(dist*dist)
                 do k=1,2
                    efield_real_1(k)=efield_real_1(k)-func_2*r_r(k)
                 enddo
                 efield_real_1(3)=efield_real_1(3)-func_2*z_r
              endif
           endif
        endif
     else
        ! removing self interaction in Madelung potential and field.
        ! if shell is not directly on core, mutual contributions are subtracted
        ! from both potential and field, otherwise only from potential.
        if(shell_model) then
           n=point_shel(j)
           if(n.ne.0) then
              if(calculate_shell) then
                 do k=1,2
                    f_r(k)=g_shel(k,n)-g_core(k,j)
                 enddo
                 z_r=f_shel(3,n)-f_core(3,j)
                 call transform_2(vector,f_r,r_r)
                 dist=sqrt(r_r(1)*r_r(1)+r_r(2)*r_r(2)+z_r*z_r)
                 if(dist.gt.1.d-7) then
                    ew_dist=ewald_const*dist
                    factor=charge_core(j)*ewald_const    
                    potential_real_0= &
                         -charge_shel(n)*two_const+ &
                         factor*fconerfc(ew_dist)
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                    factor=factor*ewald_const/dist
                    do k=1,2
                       efield_real_1(k)=efield_real_1(k)+ &
                            factor*fconerfc_1(ew_dist)*r_r(k)
                    enddo
                    efield_real_1(3)=efield_real_1(3)+ &
                         factor*fconerfc_1(ew_dist)*z_r
                 else
                    potential_real_0= &
                         -(charge_core(j)+charge_shel(n))*two_const
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                 endif
              else
                 do k=1,2
                    f_r(k)=g_core(k,j)-g_shel(k,n)
                 enddo
                 z_r=f_core(3,j)-f_shel(3,n)
                 call transform_2(vector,f_r,r_r)
                 dist=sqrt(r_r(1)*r_r(1)+r_r(2)*r_r(2)+z_r*z_r)
                 if(dist.gt.1.d-7) then
                    ew_dist=ewald_const*dist
                    factor=charge_shel(n)*ewald_const
                    potential_real_0= &
                           -charge_core(j)*two_const+ &
                           factor*fconerfc(ew_dist)
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                    factor=factor*ewald_const/dist
                    do k=1,2
                       efield_real_1(k)=efield_real_1(k)+ &
                            factor*fconerfc_1(ew_dist)*r_r(k)
                    enddo
                    efield_real_1(3)=efield_real_1(3)+ &
                         factor*fconerfc_1(ew_dist)*z_r
                 else
                    potential_real_0= &
                         -(charge_core(j)+charge_shel(n))*two_const
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                 endif
              endif
           else
              potential_real_0= &
                   -charge_core(j)*two_const
              potential_real_1=potential_real_1+ &
                   potential_real_0
           endif
        else
           potential_real_0= &
                -charge_core(j)*two_const
           potential_real_1=potential_real_1+ &
                potential_real_0
        endif
     endif
  enddo
  ! Contributions of adjacent cells in real space
  do j=1,no_core
     do ix=-max_r_x,max_r_x
        do iy=-max_r_y,max_r_y
           if(ix.ne.0.or.iy.ne.0) then
              f(1)=g_core(1,j)+ix
              f(2)=g_core(2,j)+iy
              z=f_core(3,j)
              do k=1,2
                 f_r(k)=f_probe(k)-f(k)
              enddo
              z_r=z_probe-z
              call transform_2(vector,f_r,r_r)
              dist=distance_2(vector,f_probe,z_probe,f,z)
              ! only do summation in the sphere of r_ewald_cut radius
              if(dist.le.r_ewald_cut) then
                 ew_dist=ewald_const*dist
                 func_1=dconerfc(ew_dist)/dist
                 potential_real_2=potential_real_2+ &
                      charge_core(j)*func_1
                 func_2=(two_const*exp(-ew_dist*ew_dist)+func_1)* &
                      charge_core(j)/(dist*dist)
                 do k=1,2
                    efield_real_2(k)= &
                         efield_real_2(k)-func_2*r_r(k)
                 enddo
                 efield_real_2(3)= &
                      efield_real_2(3)-func_2*z_r
                 if(shell_model) then
                    n=point_shel(j)
                    if(n.ne.0) then
                       f(1)=g_shel(1,n)+ix
                       f(2)=g_shel(2,n)+iy
                       z=f_shel(3,n)
                       do k=1,2
                          f_r(k)=f_probe(k)-f(k)
                       enddo
                       z_r=z_probe-z
                       call transform_2(vector,f_r,r_r)
                       dist=distance_2(vector,f_probe,z_probe,f,z)
                       ew_dist=ewald_const*dist
                       func_1=dconerfc(ew_dist)/dist
                       potential_real_2=potential_real_2+ &
                            charge_shel(n)*func_1
                       func_2=(two_const*exp(-ew_dist*ew_dist)+ &
                            func_1)*charge_shel(n)/(dist*dist)
                       do k=1,2
                          efield_real_2(k)=efield_real_2(k)- &
                               func_2*r_r(k)
                       enddo
                       efield_real_2(3)=efield_real_2(3)- &
                            func_2*z_r
                    endif
                 endif
              endif
           endif
        enddo
     enddo
  enddo
  ! Contribution from reciprocal lattice q vectors
  ! case y/=0
  do j=1,no_core
     do ix=-max_q_x,max_q_x
        do iy=1,max_q_y
           f_q(1)=float(ix)
           f_q(2)=float(iy)
           do k=1,2
              f_r(k)=f_probe(k)-g_core(k,j)
           enddo
           z_r=z_probe-f_core(3,j)
           call transform_2(vector,f_r,r_r)
           call transform_2(vecinv,f_q,r_q)
           qr_c=0.d0
           do k=1,2
              qr_c=qr_c+r_q(k)*r_r(k)
           enddo
           r_q_abs=sqrt(r_q(1)*r_q(1)+r_q(2)*r_q(2))
           z_abs=abs(z_r)
           arg1=r_q_abs*z_abs
           arg2=r_q_abs*0.5d0/ewald_const
           arg3=ewald_const*z_abs
           func_3 = exp(-arg1)*dconerfc(arg2-arg3)
           func_4 = -func_3
           ! Threshold avoids NaN caused by large arg1 values
           ! xbig is the value at which erfc=0 using calerf
           if (arg2 + arg3 < xbig) then
              ftmp = exp(arg1)*dconerfc(arg2+arg3)
              func_3 = func_3 + ftmp
              func_4 = func_4 + ftmp
           end if
           func_4 = func_4 * r_q_abs
           contr_c=charge_core(j)*func_3/r_q_abs
           potential_reciprocal=potential_reciprocal+ &
                contr_c*cos(qr_c)
           do k=1,2
              efield_reciprocal(k)=efield_reciprocal(k)- &
                   contr_c*sin(qr_c)*r_q(k)
           enddo
           efield_reciprocal(3)=efield_reciprocal(3)+ &
                charge_core(j)*func_4*cos(qr_c)
           if(shell_model) then
              n=point_shel(j)
              if(n.ne.0) then
                 do k=1,2
                    f_r(k)=f_probe(k)-g_shel(k,n)
                 enddo
                 call transform_2(vector,f_r,r_r)
                 qr_s=0.d0
                 do k=1,2
                    qr_s=qr_s+r_q(k)*r_r(k)
                 enddo
                 contr_s=charge_shel(n)*func_3/r_q_abs
                 potential_reciprocal=potential_reciprocal+ &
                      contr_s*cos(qr_s)
                 do k=1,2
                    efield_reciprocal(k)= &
                         efield_reciprocal(k)- &
                         contr_s*sin(qr_s)*r_q(k)
                 enddo
                 efield_reciprocal(3)=efield_reciprocal(3)+ &
                      charge_shel(n)*func_4*cos(qr_s)
              endif
           endif
           
           if (debug > 1) then 
              if(shell_model.and.n.ne.0) then
                 if(ix.eq.1.and.iy.eq.1) then
                    write(*,'(a,2g16.8/a,2g16.8/a,i2,3g20.12)') &
                         'f_q:',(f_q(k),k=1,2), &
                         'r_q:',(r_q(k),k=1,2), &
                         'j, qr_c, qr_s, r_q_2:', &
                         j,qr_c,qr_s,r_q_2, &
                         'contr_c:',contr_c
                 endif
              endif
           endif
           
        enddo
     enddo
     ! case y=0
     do ix=1,max_q_x
        f_q(1)=float(ix)
        f_q(2)=float(iy)
        do k=1,2
           f_r(k)=f_probe(k)-g_core(k,j)
        enddo
        z_r=z_probe-f_core(3,j)
        call transform_2(vector,f_r,r_r)
        call transform_2(vecinv,f_q,r_q)
        qr_c=0.d0
        do k=1,2
           qr_c=qr_c+r_q(k)*r_r(k)
        enddo
        r_q_abs=sqrt(r_q(1)*r_q(1)+r_q(2)*r_q(2))
        z_abs=abs(z_r)
        arg1=r_q_abs*z_abs
        arg2=r_q_abs*0.5d0/ewald_const
        arg3=ewald_const*z_abs
        func_3 = exp(-arg1)*dconerfc(arg2-arg3)
        func_4 = -func_3
        ! Threshold avoids NaN caused by large arg1 values
        ! xbig is the value at which erfc=0 using calerf
        if (arg2 + arg3 < xbig) then
           ftmp = exp(arg1)*dconerfc(arg2+arg3)
           func_3 = func_3 + ftmp
           func_4 = func_4 + ftmp
        end if
        func_4 = func_4 * r_q_abs
        contr_c=charge_core(j)*func_3/r_q_abs
        potential_reciprocal=potential_reciprocal+ &
             contr_c*cos(qr_c)
        do k=1,2
           efield_reciprocal(k)=efield_reciprocal(k)- &
                contr_c*sin(qr_c)*r_q(k)
        enddo
        efield_reciprocal(3)=efield_reciprocal(3)+ &
             charge_core(j)*func_4*cos(qr_c)
        if(shell_model) then
           n=point_shel(j)
           if(n.ne.0) then
              do k=1,2
                 f_r(k)=f_probe(k)-g_shel(k,n)
              enddo
              call transform_2(vector,f_r,r_r)
              qr_s=0.d0
              do k=1,2
                 qr_s=qr_s+r_q(k)*r_r(k)
              enddo
              contr_s=charge_shel(n)*func_3/r_q_abs
              potential_reciprocal=potential_reciprocal+ &
                   contr_s*cos(qr_s)
              do k=1,2
                 efield_reciprocal(k)= &
                      efield_reciprocal(k)- &
                      contr_s*sin(qr_s)*r_q(k)
              enddo
              efield_reciprocal(3)=efield_reciprocal(3)+ &
                   charge_shel(n)*func_4*cos(qr_s)
           endif
        endif
     enddo
  enddo
  ! rec_const=2*Pi/Area
  potential_reciprocal=potential_reciprocal*rec_const
  do k=1,3
     efield_reciprocal(k)=efield_reciprocal(k)*rec_const
  enddo
  potential=potential_real_1+ &
       potential_real_2+ &
       potential_reciprocal
  do k=1,3
     efield(k)=efield_real_1(k)+ &
          efield_real_2(k)+ &
          efield_reciprocal(k)
  enddo

  if (debug > 0) then  
     if (calculate_shell) then
        write(*,'(/a,i5)') "Shell of core ", j_site
     else
        write(*,'(/a,i5)') "Core ", j_site
     endif
     write(*,'(a,f26.18)') &
          'Potential (V):               ', &
          potential*ha_to_ev, &
          '   real 0 (V):               ', &
          potential_real_0*ha_to_ev, &
          '   real 1 (V):               ', &
          potential_real_1*ha_to_ev, &
          '   real 2 (V):               ', &
          potential_real_2*ha_to_ev, &
          '   reciprocal (V):           ', &
          potential_reciprocal*ha_to_ev
     write(*,'(a,3f26.18)') &
          'Derivative (V/A):            ', &
          (efield(k)*force_convert,k=1,3), &
          '   real 1  (V/A):            ', &
          (efield_real_1(k)*force_convert,k=1,3), &
          '   real 2  (V/A):            ', &
          (efield_real_2(k)*force_convert,k=1,3), &
          '   reciprocal (V/A):         ', &
          (efield_reciprocal(k)*force_convert,k=1,3)
  endif

  return
end subroutine periodic_pot_2d


subroutine periodic_pot_3d(f_probe, j_site, calculate_shell, &
  shell_model, no_core, no_shel, point_shel, &
  f_core, f_shel, charge_core, charge_shel, &
  a_cell, b_cell, c_cell, alpha_cell, beta_cell, gamma_cell, &
  vector, vecinv, volume, &
  potential, efield)
  !
  ! If j_site = 0, calculate electrostatic potential and field at point f_probe 
  ! if j_site != 0, potential and field are calculated at lattice site j_site
  !
  ! Adapted from the routine pot_3d in Construct v84
  ! 
  implicit none

  external dconerfc

  integer :: accuracy
  parameter (accuracy=20)
  double precision :: float_nr_nd
  parameter (float_nr_nd=1.1d0)

  logical :: shell_model, calculate_shell       
  double precision, dimension(3) :: f_probe(3)
  double precision :: potential, potential_real_0, potential_real_1, &
       potential_real_2, potential_reciprocal
  double precision, dimension(3) :: efield, efield_real_1, &
       efield_real_2, efield_reciprocal
  double precision, dimension(3, no_core) :: f_core, g_core
  double precision, dimension(3, no_shel) :: f_shel, g_shel
  double precision, dimension(no_core) :: charge_core
  double precision, dimension(no_shel) :: charge_shel
  
  double precision :: a_cell, b_cell, c_cell, &
       alpha_cell, beta_cell, gamma_cell
  double precision, dimension(3, 3) :: vector, vecinv
  double precision :: volume

  double precision :: pi

  double precision :: ewald_const, ewald_const_2, accuracy_factor, &
       q_const, rec_const, two_const, r_ewald_cut, r_v
  double precision :: dconerfc, fconerfc, fconerfc_1
  double precision :: distance_3, dist, ew_dist
  double precision, dimension(3) :: f, r, f_q, r_q, f_r, r_r
  double precision :: r_q_2, qr_c, qr_s, func_1, func_2, func_3
  double precision :: contr_c, contr_s, factor
  integer :: j_site, no_core, no_shel
  integer, dimension(no_core) :: point_shel
  integer :: max_r_x, max_r_y, max_r_z
  integer :: max_q_x, max_q_y, max_q_z
  integer :: ch_stdout, ix, iy, iz
  integer :: j, k, n
 
  data ch_stdout/6/

  ! Choice of extents of the supercell in reciprocal space.
  ! Seems to be independent of direct latice with aproproate choice of Ewald's
  ! constant. The current choice provides sufficient precision for real*8.
  data max_q_x/8/, max_q_y/8/, max_q_z/8/
  
  ! Initialise
  pi = dacos(-1.d0)
  potential=0.d0
  potential_real_1=0.d0
  potential_real_2=0.d0
  potential_reciprocal=0.d0
  do k=1,3
     efield(k)=0.d0
     efield_real_1(k)=0.d0
     efield_real_2(k)=0.d0
     efield_reciprocal(k)=0.d0
  enddo
  ! Choose Ewald constant and set other constants
  ewald_const_2=float_nr_nd*pi/volume**(2.d0/3.d0)
  ewald_const=sqrt(ewald_const_2)
  q_const=-1.d0/(4.d0*ewald_const_2)
  !
  ! Factor for reciprocal term is increased by a factor of 2 
  ! on collapsing complex exponents to cosine
  !
  rec_const=8.d0*pi/volume
  two_const=2.d0*ewald_const/sqrt(pi)
  !
  ! Calculate extensions of UC, which includes sphere of radius 
  ! sqrt(-ln(10^accuracy))/ewald_const
  ! 
  accuracy_factor=10.d0**accuracy
  accuracy_factor=sqrt(log(accuracy_factor))
  r_ewald_cut=accuracy_factor/ewald_const
  r_v=r_ewald_cut/volume
  max_r_x=r_v*b_cell*c_cell*abs(sin(alpha_cell))+1
  max_r_y=r_v*a_cell*c_cell*abs(sin( beta_cell))+1
  max_r_z=r_v*a_cell*b_cell*abs(sin(gamma_cell))+1
  ! IFDEBUG
  !if(j_site.eq.1) then
  !   write(*,*) 'VOLUME: ', volume
  !   write(*,'(a,g16.8,a,g16.8,a,3i4)') &
  !        'Ewald constant: ',ewald_const_2*3.57106767, &
  !        '  Ewald cut-off: ',r_ewald_cut*0.529177d0, &
  !        ' max: ',max_r_x,max_r_y,max_r_z 
     !          'q exp constant:       ',q_const, &
     !          'q sum factor:         ',rec_const
  !endif
  ! ENDDEBUG
  if(j_site.ne.0) then
     n=point_shel(j_site)
     do k=1,3
        if(calculate_shell) then
           f_probe(k)=f_shel(k,n)
        else
           f_probe(k)=f_core(k,j_site)
        endif
     enddo
  endif
  ! Centre cell around f_probe
  ! with fractional coordinates named g_core & g_shel
  do j=1,no_core
     do k=1,3
        f(k)=f_core(k,j)
     enddo
     call translate_3(vector,f_probe,f,r)
     do k=1,3
        g_core(k,j)=f(k)
        f_r(k)=f(k)-f_probe(k)
     enddo
     if(shell_model) then
        n=point_shel(j)
        if(n.ne.0) then
           do k=1,3
              g_shel(k,n)=f_shel(k,n)+(g_core(k,j)-f_core(k,j))
           enddo
        endif
     endif
  enddo
  ! 1st uc contribution in real space
  do j=1,no_core
     if(j.ne.j_site) then
        do k=1,3
           f(k)=g_core(k,j)
           f_r(k)=f_probe(k)-f(k)
        enddo
        call transform_3(vector,f_r,r_r)
        dist=distance_3(vector,f_probe,f)
        ! this check is relevant for very large UC
        if(dist.lt.r_ewald_cut) then
           ew_dist=ewald_const*dist
! IFDEBUG
!           if(j_site.eq.1) then
!              write(*,'(3x,i4,a,5g20.12)') &
!                   j,' c ',charge_core(j),dist*0.529177d0, &
!                   (f_r(k),k=1,3)
!           endif
! ENDDEBUG
           func_1=dconerfc(ew_dist)/dist
           potential_real_1=potential_real_1 + &
                charge_core(j)*func_1
           func_2=(two_const*exp(-ew_dist*ew_dist)+func_1) * &
                charge_core(j)/(dist*dist)
           do k=1,3
              efield_real_1(k)=efield_real_1(k)-func_2*r_r(k)
           enddo
           if(shell_model) then
              n=point_shel(j)
              if(n.ne.0) then
                 do k=1,3
                    f(k)=g_shel(k,n)
                    f_r(k)=f_probe(k)-f(k)
                 enddo
                 call transform_3(vector,f_r,r_r)
                 dist=distance_3(vector,f_probe,f)
! IFDEBUG
!                 write(*,'(3x,i4,a,5f16.10)') &
!                      j,' s ',charge_shel(n),dist*0.529177d0, &
!                      (f_r(k),k=1,3)
! ENDDEBUG
                 ew_dist=ewald_const*dist
                 func_1=dconerfc(ew_dist)/dist
                 potential_real_1=potential_real_1 + &
                      charge_shel(n)*func_1
                 func_2=(two_const*exp(-ew_dist*ew_dist)+func_1)* &
                      charge_shel(n)/(dist*dist)
                 do k=1,3
                    efield_real_1(k)=efield_real_1(k)-func_2*r_r(k)
                 enddo
              endif
           endif
        endif
     else
! removing self interaction in Madelung potential and field.
! if shell is not directly on core, mutual contributions are subtracted
! from both potential and field, otherwise only from potential.
        if(shell_model) then
           n=point_shel(j)
           if(n.ne.0) then
              if(calculate_shell) then
                 do k=1,3
                    f_r(k)=g_shel(k,n)-g_core(k,j)
                 enddo
                 call transform_3(vector,f_r,r_r)
                 dist=sqrt(r_r(1)*r_r(1)+r_r(2)*r_r(2)+r_r(3)*r_r(3))
                 if(dist.gt.1.d-7) then
                    ew_dist=ewald_const*dist
                    factor=charge_core(j)*ewald_const                    
                    potential_real_0= &
                         -charge_shel(n)*two_const+ &
                         factor*fconerfc(ew_dist)
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                    factor=factor*ewald_const/dist
                    do k=1,3
                       efield_real_1(k)=efield_real_1(k)+ &
                            factor*fconerfc_1(ew_dist)*r_r(k)
                    enddo
                 else
                    potential_real_0= &
                         -(charge_core(j)+charge_shel(n))*two_const
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                 endif
              else
                 do k=1,3
                    f_r(k)=g_core(k,j)-g_shel(k,n)
                 enddo
                 call transform_3(vector,f_r,r_r)
                 dist=sqrt(&
                      r_r(1)*r_r(1)+r_r(2)*r_r(2)+r_r(3)*r_r(3))
                 if(dist.gt.1.d-7) then
                    ew_dist=ewald_const*dist
                    factor=charge_shel(n)*ewald_const
                    potential_real_0= &
                         -charge_core(j)*two_const+ &
                         factor*fconerfc(ew_dist)
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                    factor=factor*ewald_const/dist
                    do k=1,3
                       efield_real_1(k)=efield_real_1(k)+ &
                            factor*fconerfc_1(ew_dist)*r_r(k)
                    enddo
                 else
                    potential_real_0= &
                         -(charge_core(j)+charge_shel(n))*two_const
                    potential_real_1=potential_real_1+ &
                         potential_real_0
                 endif
              endif
           else
              potential_real_0= &
                   -charge_core(j)*two_const
              potential_real_1=potential_real_1+ &
                   potential_real_0
           endif
        else
           potential_real_0= &
                -charge_core(j)*two_const
           potential_real_1=potential_real_1+ &
                potential_real_0
        endif
     endif
  enddo
! Contributions of adjacent cells in real space
  do j=1,no_core
     do ix=-max_r_x,max_r_x
        do iy=-max_r_y,max_r_y
           do iz=-max_r_z,max_r_z
              if(ix.ne.0.or.iy.ne.0.or.iz.ne.0) then
                 f(1)=g_core(1,j)+ix
                 f(2)=g_core(2,j)+iy
                 f(3)=g_core(3,j)+iz
                 do k=1,3
                    f_r(k)=f_probe(k)-f(k)
                 enddo
                 call transform_3(vector,f_r,r_r)
                 dist=distance_3(vector,f_probe,f)
! only do summation in the sphere of r_ewald_cut radius
                 if(dist.le.r_ewald_cut) then
                    ew_dist=ewald_const*dist
                    func_1=dconerfc(ew_dist)/dist
                    potential_real_2=potential_real_2+ &
                         charge_core(j)*func_1
                    func_2=(two_const*exp(-ew_dist*ew_dist)+func_1)* &
                         charge_core(j)/(dist*dist)
                    do k=1,3
                       efield_real_2(k)= &
                            efield_real_2(k)-func_2*r_r(k)
                    enddo
                    if(shell_model) then
                       n=point_shel(j)
                       if(n.ne.0) then
                          f(1)=g_shel(1,n)+ix
                          f(2)=g_shel(2,n)+iy
                          f(3)=g_shel(3,n)+iz
                          do k=1,3
                             f_r(k)=f_probe(k)-f(k)
                          enddo
                          call transform_3(vector,f_r,r_r)
                          dist=distance_3(vector,f_probe,f)
                          ew_dist=ewald_const*dist
                          func_1=dconerfc(ew_dist)/dist
                          potential_real_2=potential_real_2+ &
                               charge_shel(n)*func_1
                          func_2=(two_const*exp(-ew_dist*ew_dist)+ &
                               func_1)*charge_shel(n)/(dist*dist)
                          do k=1,3
                             efield_real_2(k)=efield_real_2(k)- &
                                  func_2*r_r(k)
                          enddo
                       endif
                    endif
                 endif
              endif
           enddo
        enddo
     enddo
  enddo
! Contribution from reciprocal lattice q vectors
! with zero components along x and y axes and z=0, treated separately
  do j=1,no_core
     do ix=-max_q_x,max_q_x
        do iy=-max_q_y,max_q_y
           do iz=1,max_q_z
              f_q(1)=float(ix)
              f_q(2)=float(iy)
              f_q(3)=float(iz)
              do k=1,3
                 f_r(k)=f_probe(k)-g_core(k,j)
              enddo
              call transform_3(vector,f_r,r_r)
              call transform_3(vecinv,f_q,r_q)
              qr_c=0.d0
              do k=1,3
                 qr_c=qr_c+r_q(k)*r_r(k)
              enddo
              r_q_2=r_q(1)*r_q(1)+r_q(2)*r_q(2)+r_q(3)*r_q(3)
              func_3=exp(q_const*r_q_2)/r_q_2
              contr_c=charge_core(j)*func_3
              potential_reciprocal=potential_reciprocal+ &
                   contr_c*cos(qr_c)
              do k=1,3
                 efield_reciprocal(k)=efield_reciprocal(k)- &
                      contr_c*sin(qr_c)*r_q(k)
              enddo
              if(shell_model) then
                 n=point_shel(j)
                 if(n.ne.0) then
                    do k=1,3
                       f_r(k)=f_probe(k)-g_shel(k,n)
                    enddo
                    call transform_3(vector,f_r,r_r)
                    qr_s=0.d0
                    do k=1,3
                       qr_s=qr_s+r_q(k)*r_r(k)
                    enddo
                    contr_s=charge_shel(n)*func_3
! q_const=-1/(2*Ewald_const)**2
                    potential_reciprocal=potential_reciprocal+ &
                         contr_s*cos(qr_s)
                    do k=1,3
                       efield_reciprocal(k)= &
                            efield_reciprocal(k)- &
                            contr_s*sin(qr_s)*r_q(k)
                    enddo
                 endif
              endif
! IFDEBUG
!                  if(shell_model.and.n.ne.0) then
!                     if(ix.eq.0.and.iy.eq.0.and.iz.eq.1) then
!                        write(*,'(a,3g16.8/a,3g16.8/a,i2,3g20.12)')
!     &                       'f_q:',(f_q(k),k=1,3),
!     &                       'r_q:',(r_q(k),k=1,3),
!     &                       'j, qr_c, qr_s, r_q_2:',
!     &                       j,qr_c,qr_s,r_q_2,
!     &                       'contr_c,contr_s:',contr_c,contr_s
!                     endif
!                  endif
! ENDDEBUG
           enddo
        enddo
     enddo
! z=0 case y/=0
     do ix=-max_q_x,max_q_x
        do iy=1,max_q_y
           f_q(1)=float(ix)
           f_q(2)=float(iy)
           f_q(3)=0.d0
           do k=1,3
              f_r(k)=f_probe(k)-g_core(k,j)
           enddo
           call transform_3(vector,f_r,r_r)
           call transform_3(vecinv,f_q,r_q)
           qr_c=0.d0
           do k=1,3
              qr_c=qr_c+r_q(k)*r_r(k)
           enddo
           r_q_2=r_q(1)*r_q(1)+r_q(2)*r_q(2)+r_q(3)*r_q(3)
           func_3=exp(q_const*r_q_2)/r_q_2
           contr_c=charge_core(j)*func_3
           potential_reciprocal=potential_reciprocal+ &
                contr_c*cos(qr_c)
           do k=1,3
              efield_reciprocal(k)=efield_reciprocal(k)- &
                   contr_c*sin(qr_c)*r_q(k)
           enddo
           if(shell_model) then
              n=point_shel(j)
              if(n.ne.0) then
                 do k=1,3
                    f_r(k)=f_probe(k)-g_shel(k,n)
                 enddo
                 call transform_3(vector,f_r,r_r)
                 qr_s=0.d0
                 do k=1,3
                    qr_s=qr_s+r_q(k)*r_r(k)
                 enddo
                 contr_s=charge_shel(n)*func_3
! q_const=-1/(2*Ewald_const)**2
                 potential_reciprocal=potential_reciprocal+ &
                      contr_s*cos(qr_s)
                 do k=1,3
                    efield_reciprocal(k)= &
                         efield_reciprocal(k)- &
                         contr_s*sin(qr_s)*r_q(k)
                 enddo
              endif
           endif
! IFDEBUG
!               if(shell_model.and.n.ne.0) then
!                  if(ix.eq.1.and.iy.eq.1) then
!                     write(*,'(a,3g16.8/a,3g16.8/a,i2,3g20.12)')
!     &                    'f_q:',(f_q(k),k=1,3),
!     &                    'r_q:',(r_q(k),k=1,3),
!     &                    'j, qr_c, qr_s, r_q_2:',
!     &                    j,qr_c,qr_s,r_q_2,
!     &                    'contr_c:',contr_c
!                  endif
!               endif
! ENDDEBUG
        enddo
     enddo
! z=0 and y=0
     do ix=1,max_q_x
        f_q(1)=float(ix)
        f_q(2)=0.d0
        f_q(3)=0.d0
        do k=1,3
           f_r(k)=f_probe(k)-g_core(k,j)
        enddo
        call transform_3(vector,f_r,r_r)
        call transform_3(vecinv,f_q,r_q)
        qr_c=0.d0
        do k=1,3
           qr_c=qr_c+r_q(k)*r_r(k)
        enddo
        r_q_2=r_q(1)*r_q(1)+r_q(2)*r_q(2)+r_q(3)*r_q(3)
        func_3=exp(q_const*r_q_2)/r_q_2
        contr_c=charge_core(j)*func_3
        potential_reciprocal=potential_reciprocal+ &
             contr_c*cos(qr_c)
        do k=1,3
           efield_reciprocal(k)=efield_reciprocal(k)- &
                contr_c*sin(qr_c)*r_q(k)
        enddo
        if(shell_model) then
           n=point_shel(j)
           if(n.ne.0) then
              do k=1,3
                 f_r(k)=f_probe(k)-g_shel(k,n)
              enddo
              call transform_3(vector,f_r,r_r)
              qr_s=0.d0
              do k=1,3
                 qr_s=qr_s+r_q(k)*r_r(k)
              enddo
              contr_s=charge_shel(n)*func_3
 ! q_const=-1/(2*Ewald_const)**2
              potential_reciprocal=potential_reciprocal+ &
                   contr_s*cos(qr_s)
              do k=1,3
                 efield_reciprocal(k)= &
                      efield_reciprocal(k)- &
                      contr_s*sin(qr_s)*r_q(k)
              enddo
           endif
        endif
     enddo
  enddo
! rec_const=8.d0*Pi/Volume
  potential_reciprocal=potential_reciprocal*rec_const
  do k=1,3
     efield_reciprocal(k)=efield_reciprocal(k)*rec_const
  enddo
  potential=potential_real_1+ &
       potential_real_2+ &
       potential_reciprocal
  do k=1,3
     efield(k)=efield_real_1(k)+ &
          efield_real_2(k)+ &
          efield_reciprocal(k)
  enddo
  return
end subroutine periodic_pot_3d
