diff --git a/Na_Cl.in b/Na_Cl.in
new file mode 100644
index 0000000..b9faa06
--- /dev/null
+++ b/Na_Cl.in
@@ -0,0 +1,15 @@
+NaCl
+  1.0
+  1.0 1.0 1.0
+  200  200  200
+  8
+ 1.0q0  0.0q0  0.0q0  0.0q0
+ 1.0q0  0.5q0  0.0q0  0.5q0
+ 1.0q0  0.0q0  0.5q0  0.5q0
+ 1.0q0  0.5q0  0.5q0  0.0q0
+-1.0q0  0.5q0  0.0q0  0.0q0
+-1.0q0  0.0q0  0.5q0  0.0q0
+-1.0q0  0.0q0  0.0q0  0.5q0
+-1.0q0  0.5q0  0.5q0  0.5q0
+ 2.0
+END
diff --git a/madelung.f90 b/madelung.f90
new file mode 100644
index 0000000..e118a3c
--- /dev/null
+++ b/madelung.f90
@@ -0,0 +1,103 @@
+program madelung
+
+implicit none
+
+external :: dist
+integer :: nx, ny, nz, ix, iy, iz
+integer :: nat, iat, jat
+real*8 :: r0, r1, r2, rexp, fact
+real*8 :: dx, dy, dz
+real*8 :: qtot, etot, e0, e1, e2, e3
+real*8, allocatable :: q(:), x(:), y(:), z(:)
+real*8 :: xx, yy, zz
+real*8 :: dist, dd
+real*8 :: start, finish
+r0=0.0_8
+r1=1.0_8
+r2=2.0_8
+!
+!Input (free format)
+!
+!	Title
+!	rexp		(minus) potential exponent (1.0 for Coulomb)
+!	dx, dy, dz	unit-cell dimensions (in bohr)
+!	nx, ny, nz	unit-cell number
+!	nat		number of centers in the unit cell
+!	loop over nat 
+!		q x y z		atom charge (in a.u.), atom position (in bohr)
+!	end loop
+!	fact		a scale factor to divide the final energy
+!	END
+!
+read (*,*)
+read (*,*) rexp
+read (*,*) dx, dy, dz
+read (*,*) nx, ny, nz
+read (*,*) nat
+allocate (q(1:nat),x(1:nat),y(1:nat),z(1:nat))
+do iat=1,nat
+   read (5,*) q(iat), x(iat), y(iat), z(iat)
+enddo
+read (*,*) fact
+!
+call cpu_time(start)
+!
+qtot=r0
+etot=r0
+do iat=1,nat                              ! test cell
+   qtot=qtot+q(iat)
+   e0=r0
+   do ix=0,nx-1
+      e1=r0
+      do iy=0,ny-1
+         e2=r0
+         do iz=0,nz-1
+            e3=r0
+            do jat=1,nat                  ! general cell
+               xx=real(ix,8)+x(iat)-x(jat)
+               yy=real(iy,8)+y(iat)-y(jat)
+               zz=real(iz,8)+z(iat)-z(jat)
+               dd=dist(xx,yy,zz,nx,ny,nz)
+!                    write (6,*) iat, jat, ix, iy, iz, dd
+!                    if (dd < 1.q-10) cycle
+               if (ix==0.and.iy==0.and.iz==0.and.iat==jat) cycle
+               e3=e3+q(iat)*q(jat)/dd**rexp
+            enddo
+            e2=e2+e3
+         enddo
+         e1=e1+e2
+      enddo
+      e0=e0+e1
+   enddo
+   etot=etot+e0
+enddo
+!
+call cpu_time(finish)
+!
+write (6,*) 
+write (6,101) qtot             ! total charge - is a check: for neutral systems qtot should be zero
+write (6,102) etot             ! unit-cell energy / nat
+write (6,103) etot/real(nat,8) ! unit-cell energy / (nat*fact)
+write (6,104) etot/real(nat,8)/fact 
+write (6,*) 
+print '("Time = ",f6.3," seconds.")',finish-start
+stop
+101   format (f20.12,'	total charge')
+102   format (f20.12,'	unit-cell energy')
+103   format (f20.12,'	per-atom cell energy')
+104   format (f20.12,'	per-atom cell energy / fact')
+end program madelung
+
+function dist(x,y,z,nx,ny,nz)
+implicit none
+integer :: nx, ny, nz
+real*8 :: x, y, z
+real*8 :: r1, pi
+real*8 :: dist
+r1=1.0_8
+pi=acos(-r1)
+dist=((real(nx,8)/pi*sin(pi*x/real(nx,8)))**2+(real(ny,8)/pi*sin(pi*y/real(ny,8)))**2 &
+     +(real(nz,8)/pi*sin(pi*z/real(nz,8)))**2)
+dist=sqrt(dist)
+return
+end function dist