PROGRAM MADEL
************************************************************************
*                CALCULATE MADELUNG POTENTIAL
*                      V1.02  14.10.88
************************************************************************

* MADEL, VER. 1.02A
*     THIS PROGRAM WAS CODED BY K. KATO, NATIONAL INSTITUTE FOR RESEARCH
*     IN INORGANIC MATERIAL AND SLIGHTLY MODIFIED BY F. IZUMI ON
*     13 JULY, 1990.

      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (NNA=100, P=3.141592653589793D0)
      PARAMETER (NSP=200000)
*
      REAL*8 FA(NSP),FB(NSP)
      COMMON /COM1/ FA,FB
      COMMON IZ,NS,NA,RADIUS,REGION
      COMMON RCP(6),X(3,NNA),Z(NNA),W(NNA)
      COMMON TS(3,48),IS(2,3,48)
      CHARACTER ATOM(NNA)*8, TITEL*72
*SUN
      CHARACTER*50 FILE5
      REAL*8 X1(3),X2(3),PATOMP(NNA)
*
      WRITE (6,'(1H1,A/1H ,A/1H ,A/)')
     &  '***************************************',
     &  '* MADELUNG-POTENTIAL (V1.02 14.10.88) *',
     &  '***************************************'
*SUN
      CALL GETARG(1,FILE5)
*SUN
      OPEN(UNIT=5,FILE=FILE5,STATUS='OLD')
      READ (5,'(A72)') TITEL
      WRITE (6,'(1H ,A72)') TITEL
      READ (5,'(2F6.2)') RADIUS, REGION
      WRITE (6,2010) RADIUS,REGION
 2010 FORMAT (' RADIUS DER IONENKUGEL',F6.2/
     &  ' SUMMATIONSBEREICH IM REZIPROKEN RAUM (1/ANGSTROM)',F6.2)
      READ (5,'(3I3)') NS,NA,IZ
      WRITE (6,'('' ANZAHL DER SYM.OPERATIONEN'',I4)') NS
      WRITE (6,'('' ANZAHL DER ATOME'',I4)') NA
      WRITE (6,'('' SYMMETRIEZENTRUM'',I4)') IZ
*
      READ (5,'(6F9.6)') A,B,C,AL,BE,GA
      WRITE (6,'(/'' GITTERKONSTANTEN '',3F9.4,3F9.3)') A,B,C,AL,BE,GA
*
      READ (5,'(F11.6,2I2,F11.6,2I2,F11.6,2I2)')
     &    ((TS(I,J),(IS(K,I,J),K=1,2),I=1,3),J=1,NS)
      WRITE (6,'(/'' SYMMETRIEOPERATIONEN''/)')
      WRITE (6,'(1H ,F13.6,2I2,F13.6,2I2,F13.6,2I2)')
     &    ((TS(I,J),(IS(K,I,J),K=1,2),I=1,3),J=1,NS)
*
      READ (5,'(A8,1X,5F9.4)') (ATOM(I),Z(I),W(I),(X(J,I),J=1,3),I=1,NA)
*
      CALL SLTP (A,B,C,AL,BE,GA,V,RA,RB,RC,RAL,RBE,RGA)
      IHX = A * REGION
      IKX = B * REGION
      ILX = C * REGION
*
      IHKL = (2*IHX+1)*(2*IKX+1)*(ILX+1)
      WRITE (6,'(/'' SPEICHERBEDARF'',3I6,I10,'' ('',I6,'')'')')
     &  IHX,IKX,ILX,IHKL,NSP
      IF (IHKL .GT. NSP) STOP 1
*
      RCP(1) = RA * RA
      RCP(2) = RB * RB
      RCP(3) = RC * RC
      RCP(4) = RA * RB * RGA * 2.0D0
      RCP(5) = RA * RC * RBE * 2.0D0
      RCP(6) = RB * RC * RAL * 2.0D0
*
      CALL FC (FA,FB,IHX,IKX,ILX)
*
*     -- Potential an den Atomlagen
*
      WRITE (6,2020)
 2020 FORMAT (/' POTENTIAL AN DEN ATOMLAGEN'//12X,'LADUNG    W',8X,'X',
     &  8X,'Y',8X,'Z       POTENTIAL'/1H )
      DO 30 I = 1, NA
        PATOMP(I) = FOURR (FA,FB,IHX,IKX,ILX,IZ,X(1,I)) / P / V
     &           - 21.0D0 * Z(I) / RADIUS / 8.0D0
        WRITE (6,'(1H ,A8,5F9.6,1P,E15.6)')
     &     ATOM(I),Z(I),W(I),(X(J,I),J=1,3),PATOMP(I)
   30 CONTINUE
*
*     -- Elektrostatische Energie der asymmetrischen Einheit
*
      EE = 0.0D0
      DO 40 I = 1, NA
        EE = Z(I) * W(I) * PATOMP(I) + EE
   40 CONTINUE
      EE = EE * 0.5D0
      WRITE (6,2022) EE, 14.39975D0*EE, 1.389354D0*EE, 331.901D0*EE
 2022 FORMAT (//,' ELEKTROSTATISCHE ENERGIE DER ASYMMETRISCHEN EINHEIT'/
* MODIFIED BY F. IZUMI, JULY 10, 1990
*    &  3X,F11.4,' E**2/A,',F11.3,' EV,',F11.4,' MJ/MOL,',F11.1,
     &  3X,F13.6,' E**2/A,',F13.5,' EV,',F13.6,' MJ/MOL,',F13.3,
     &  ' KCAL/MOL')
*
*     -- Potential an beliebigem Punkt
*
      G1 = A * A
      G2 = B * B
      G3 = C * C
      G4 = A * B * GA * 2.0D0
      G5 = A * C * BE * 2.0D0
      G6 = B * C * AL * 2.0D0
*
*Added by F. Izumi to suppress extra output
      READ (5,'(3F9.6)',END=150) X1
      BACKSPACE 5
      WRITE (6,2030)
 2030 FORMAT (/' POTENTIAL AN BELIEBIGEN PUNKTEN'//5X,'X',8X,'Y',
     &  8X,'Z',9X,'POBBA'10X,'PMBBA',10X,'B.ATOM   ABSTAND'/1H )
*
*     -- POBBA: Potential ohne Beitrag des benachbarten Atoms
*     -- PMBBA: Potential mit  Beitrag des benachbarten Atoms
*
  100 CONTINUE
      READ (5,'(3F9.6)',END=150) X1
      PMBBA = FOURR (FA,FB,IHX,IKX,ILX,IZ,X1) / P / V
      R = RADIUS * 100.0D0
*
      DO 117 IA = 1, NA
        DO 116 JS = 1, NS
          DO 112 I = 1, 3
            X2(I) = TS(I,JS)
            DO 111 J = 1, 2
              M = IS(J,I,JS)
              IF (M .GT. 0) THEN
                X2(I) = X2(I) + X(M,IA)
              ELSE IF (M .LT. 0) THEN
                X2(I) = X2(I) - X(-M,IA)
              END IF
  111       CONTINUE
  112     CONTINUE
          DO 115 IIX = -1, 1
            U1 = X1(1) - X2(1) + IIX
            ABST1 = U1 * U1 * G1
            DO 114 IIY = -1, 1
              U2 = X1(2) - X2(2) + IIY
              ABST2 = (U2 * G2 + U1 * G4) * U2 + ABST1
              DO 113 IIZ = -1, 1
                U3 = X1(3) - X2(3) + IIZ
                ABST3 = (U3 * G3 + U1 * G5 + U2 * G6) * U3 + ABST2
                ABST = SQRT (ABST3)
                IF (ABST .LT. R) THEN
                  R = ABST
                  IAMIN = IA
                  IF (R .LT. RADIUS) GO TO 118
                END IF
  113         CONTINUE
  114       CONTINUE
  115     CONTINUE
  116   CONTINUE
  117 CONTINUE
      GO TO 120
*
*     -- Der Punkt liegt innerhalb der Kugel
*
  118 CONTINUE
      IF (R .LT. 0.00001) THEN
        WRITE (6,'(1H ,3F9.6)') X1
        GO TO 100
      ELSE
        RS = R / RADIUS
        POBBA = PMBBA - (210.0D0 * Z(IA) / RADIUS / 8.0D0) *
     &    (RS**6/14.0D0 - RS**5*4.0D0/15.0D0 + RS**4*0.3D0 - 
     &    RS**2/6.0D0 + 0.1D0)
        PMBBA = POBBA + Z(IA) / R
      END IF
      GO TO 130
*
*     -- Der Punkt liegt ausserhalb der Kugel
*
  120 CONTINUE
      POBBA = PMBBA - Z(IAMIN) / R
  130 CONTINUE
      WRITE (6,'(1H ,3F9.6,1P,2E15.6,5X,A8,0P,F8.4)')
     &  X1,POBBA,PMBBA,ATOM(IAMIN),R
      GO TO 100
*
  150 CONTINUE
      STOP
      END
*
      SUBROUTINE SLTP (A,B,C,CA,CB,CG,V,RA,RB,RC,RCA,RCB,RCG)
************************************************************************
*
*       READ REAL/RECIPROCAL LATTICE PARAMETERS AND CALCULATE
*               RECIPROCAL/REAL LATTICE PARAMETERS
*                      IN SIGLE PRECISION
*
*     ---- INPUT DATA ----
*
*     A,B,C,CA,CB,CG...REAL OR RECIPROCAL LATTICE PARAMETERS
*                      WITH ANGLES OR THEIR COSINES
*
*     ---- OUTPUT PARAMETER ---
*
*     A,B,C,CA,CB,CG...REAL LATTICE PARAMETERS (CA = COS (ALPHA) ETC)
*     V...VOLUME OF THE REAL LATTICE UNIT CELL
*     RA,RB,RC,RCA,RCB,RCG...RECIPROCAL LATTICE PARAMETERS
*                            (RCA = COS (ALPHA*) ETC)
*
************************************************************************
*
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (P=3.141592653589793D0/180.0D0)
      PARAMETER (ONE = 1.0D0)
*
      IF (CA .GE. 1.0) THEN
        CA=COS(P*CA)
        CB=COS(P*CB)
        CG=COS(P*CG)
      END IF
      V=A*B*C*SQRT(1.0D0-CA**2-CB**2-CG**2+2.0D0*CA*CB*CG)
      SA=SQRT(1.0D0-CA**2)
      SB=SQRT(1.0D0-CB**2)
      SG=SQRT(1.0D0-CG**2)
      RA=B*C*SA/V
      RB=C*A*SB/V
      RC=A*B*SG/V
      RCA=(CB*CG-CA)/SB/SG
      RCB=(CG*CA-CB)/SG/SA
      RCG=(CA*CB-CG)/SA/SB
*
      IF (A .LT. 1.0) THEN
        T = A
        A = RA
        RA = T
        T = B
        B = RB
        RB = T
        T = C
        C = RC
        RC = T
        T = CA
        CA = RCA
        RCA = T
        T = CB
        CB = RCB
        RCB = T
        T = CG
        CG = RCG
        RCG = T
        V = ONE / V
      END IF
*
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION FF (R, S, I)
************************************************************************
*                                FF
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (P=3.141592653589793D0, PP=2.0D0*P)
      ARG = PP * R * S
      GO TO (10, 20) I
   10 CONTINUE
*
*     -- RHO0*(1.0 - (R/S)**2)
*
      FF = ( (3.0D0 - ARG**2) * SIN (ARG) - 3.0D0 * ARG * COS (ARG) )
     &         * 15.0D0 / ARG**5
      GO TO 100
*
   20 CONTINUE
*
*     -- RHO0*(1.0 - 6.0*(R/S)**2  + 8.0*(R/S)**3 - 3.0*(R/S)**4))
*
      FF = ((ARG**2-15.0D0)*SIN(ARG) + 7.0D0*ARG*COS(ARG) + 8.0D0*ARG)
     &         * 630.0D0 / ARG**7
*
  100 CONTINUE
      RETURN
      END
      SUBROUTINE FC (FA, FB, IHX, IKX, ILX)
***********************************************************************
*                                FC
***********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (NNA=100, PP=3.141592653589793D0*2.0D0)
      COMMON IZ,NS,NA,RADIUS,REGION
      COMMON RCP(6),X(3,NNA),Z(NNA),W(NNA)
      COMMON TS(3,48),IS(2,3,48)
      REAL*8 FA(0-IHX:IHX,0-IKX:IKX,0:ILX)
      REAL*8 FB(0-IHX:IHX,0-IKX:IKX,0:ILX)
      REAL*8 H(3),HJ(3)
*
      IKN = 0
      IHN = 1
      DO 230 IL = 0, ILX
        H(3) = IL
        DO 220 IK = IKN, IKX
          H(2) = IK
          DO 210 IH = IHN, IHX
            H(1) = IH
*
*     -- Strukturfaktor
*
      Q = H(1)**2*RCP(1) + H(2)**2*RCP(2) + H(3)**2*RCP(3)
     &      + H(1)*H(2)*RCP(4) + H(1)*H(3)*RCP(5) + H(2)*H(3)*RCP(6)
      SQRTQ = SQRT (Q)
      IF (SQRTQ .GT. REGION) THEN
        FA(IH,IK,IL) = 0.0D0
        FB(IH,IK,IL) = 0.0D0
        GO TO 210
      END IF
      FFI = FF (SQRTQ, RADIUS, 2)
*
      A = 0.0D0
      B = 0.0D0
*
      DO 150 J = 1, NS
        TJ = 0.0D0
        DO 36 I = 1, 3
          TJ = TJ + H(I) * TS(I,J)
          HJ(I) = 0.0D0
          DO 35 K = 1, 3
            DO 34 L = 1, 2
              M = IS(L,K,J)
              IF (M-I) 32,31,32
   31         HJ(I) = HJ(I) + H(K)
              GO TO 34
   32         IF (M+I) 34,33,34
   33         HJ(I) = HJ(I) - H(K)
   34       CONTINUE
   35     CONTINUE
   36   CONTINUE
*
        DO 120 I = 1, NA
          HX = TJ
          DO 110 K = 1, 3
            HX = HX + HJ(K) * X(K,I)
  110     CONTINUE
          HX = HX * PP
          A =  Z(I) * FFI * W(I) * COS (HX) + A
          IF (IZ .EQ. 0) B = Z(I) * FFI * W(I) * SIN (HX) + B
  120   CONTINUE
*
  150 CONTINUE
*
*     -- Koeffizienten der Fourierreihe des Potentials
*
      IF (IZ .EQ. 0) THEN
        FA(IH,IK,IL) = A / Q
        FB(IH,IK,IL) = B / Q
      ELSE
        FA(IH,IK,IL) = 2.0D0 * A / Q
        FB(IH,IK,IL) = 0.0D0
      END IF
*
*
*
  210     CONTINUE
          IHN = -IHX
  220   CONTINUE
        IKN = -IKX
  230 CONTINUE
*
      RETURN
      END

      DOUBLE PRECISION FUNCTION FOURR (FA, FB, IHX, IKX, ILX, IZ, X)
************************************************************************
*                        FOURIERSYNTHESE
************************************************************************
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (PP=3.141592653589793D0*2.0D0)
      REAL*8 FA(0-IHX:IHX,0-IKX:IKX,0:ILX)
      REAL*8 FB(0-IHX:IHX,0-IKX:IKX,0:ILX)
      REAL*8 X(3)
*
      X1 = X(1) * PP
      X2 = X(2) * PP
      X3 = X(3) * PP
*
      FF = 0.0D0
      IKN = 0
      IHN = 1
      DO 33 IL = 0, ILX
        DO 32 IK = IKN, IKX
          DO 31 IH = IHN, IHX
            FF = FA(IH,IK,IL) * COS (IL * X3 + IK * X2 + IH * X1) + FF
   31     CONTINUE
          IHN = -IHX
   32   CONTINUE
        IKN = -IKX
   33 CONTINUE
*
      IF (IZ .EQ. 0) THEN
        IKN = 0
        IHN = 1
        DO 36 IL = 0, ILX
          DO 35 IK = IKN, IKX
            DO 34 IH = IHN, IHX
              FF = FB(IH,IK,IL) * SIN (IL * X3 + IK * X2 + IH * X1) + FF
   34       CONTINUE
            IHN = -IHX
   35     CONTINUE
          IKN = -IKX
   36   CONTINUE
      END IF
*
      FOURR = FF * 2.0D0
*
      RETURN
      END