C====================================================
      SUBROUTINE ESANI (ESA, WHICH)
C
C Calls ESANI2, which does the actual energy calculation
C
C by Nico Tjandra and Marius Clore June 1997
C====================================================
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'susc_anis.inc'
      INCLUDE 'heap.inc'
C i/o
      DOUBLE PRECISION ESA
      CHARACTER*7 WHICH
C
      IF (SANIIPTR.NE.0) THEN
      CALL ESANI2(ESA, HEAP(SANIIPTR), HEAP(SANIJPTR), HEAP(SANIKPTR),
     &        HEAP(SANILPTR), HEAP(SANIMPTR), HEAP(SANINPTR),
     &        HEAP(SANIOBSPTR), HEAP(SANIERRPTR),
     &        HEAP(CALCSANIPTR), WHICH)
      END IF
      RETURN
      END
C
C=====================================================
      SUBROUTINE ESANI2 (ESA, ATOMI, ATOMJ, ATOMK, ATOML,
     & ATOMM, ATOMN, SANIOBS, SANIERR, SANICALC, WHICH)
C
C Calculates anisotropy energies
C
C Anis energies are of the form
C      E = k1*deltaAnis**2
C where
C      k1 = force constant,
C      deltaAnis = calculated Anis - observed Anis
C Alos note: that one defines the Anisotropy functions as:
C dip=coef0+coef1*(3*cos(theta)^2-1)+(3/2)*coef2*(sin(theta)^2*cos(2*phi))
C
C which is a flag that switches between energy & force and
C calculated Anis (for violations) calcs
C
C by Nico Tjandra and Marius Clore June 1997
C
C=========================================================
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'numbers.inc'
      INCLUDE 'deriv.inc'
      INCLUDE 'susc_anis.inc'
      INCLUDE 'consta.inc'
      INCLUDE 'mtf.inc'
C i/o
      INTEGER ATOMI(*), ATOMJ(*), ATOMK(*), ATOML(*)
      INTEGER ATOMM(*), ATOMN(*)
      DOUBLE PRECISION SANIOBS(*), SANIERR(*), SANICALC(*)
      DOUBLE PRECISION ESA
      CHARACTER*7 WHICH
C local variables
      INTEGER COUNT, CLASS
      DOUBLE PRECISION XI, XJ, XK, XL, XM, XN, YI, YJ, YK, YL, YM, YN,
     &                 ZI, ZJ, ZK, ZL, ZM, ZN,
     &                 CALCSANI, DELTASANI, E, DF1, DF2,
     &                 PHI, DPXI, DPYI, DPZI, DPXJ, DPYJ, DPZJ,
     &                 DPXK, DPYK, DPZK, DPXL, DPYL, DPZL,
     &                 DPXM, DPYM, DPZM, DPXN, DPYN, DPZN,
     &                 THETA, DTXI, DTYI, DTZI, DTXJ, DTYJ, DTZJ,
     &                 DTXK, DTYK, DTZK, DTXL, DTYL, DTZL,
     &                 DTXM, DTYM, DTZM, DTXN, DTYN, DTZN,
     &                 o1, o2, o3, o4 ,o5 ,o6 ,o7 ,o8, o9, o10,
     &                 o11, o12, o13, o14, o15, o16, o17, o18, o19,
     &                 o20, o21, o22, o23, o24, o25, o26, o27, o28,
     &                 o29, o30, o31, o32, o33, o34, o35, o36, o37,
     &                 o38, o39, o40, o41, o42, o43, o44, o45, o46,
     &                 o47, o48, o49, o50, o51, o52, o53, o54, o55,
     &                 o56, o57, o58, o59, o60, o61, o62, o63, o64,
     &                 o65, o66, o67, o68, o69, o70, o71, o72, o73,
     &                 o74, o75, o76, o77, o78, o79, o80, o81, o82
      DOUBLE PRECISION OBSSANI, ERRSANI, K1, COEF0,
     &                 COEF1, COEF2, A, B
C
C begin
C
C following Axel's code in ETOR,
C
C zero out partial energy
C
      ESA = ZERO
C
      CLASS = 1
      K1 = SANIFORCES(1)
      COEF0 = SANICOEF0(1)
      COEF1 = SANICOEF1(1)
      COEF2 = SANICOEF2(1)
C
          COUNT=0
C
      DO WHILE(COUNT.LT.NSANI)
          COUNT = COUNT+1
C
          IF (SANIASSNDX(CLASS).LT.COUNT) THEN
               CLASS = CLASS + 1
             IF (SANIASSNDX(CLASS).EQ.SANIASSNDX(CLASS-1)) THEN
               COUNT = COUNT - 1
             ELSE
               K1 = SANIFORCES(CLASS)
               COEF0 = SANICOEF0(CLASS)
               COEF1 = SANICOEF1(CLASS)
               COEF2 = SANICOEF2(CLASS)
            END IF
          END IF     
C
          XI = X(ATOMI(COUNT))
          XJ = X(ATOMJ(COUNT))
          XK = X(ATOMK(COUNT))
          XL = X(ATOML(COUNT))
          XM = X(ATOMM(COUNT))
          XN = X(ATOMN(COUNT))
          YI = Y(ATOMI(COUNT))
          YJ = Y(ATOMJ(COUNT))
          YK = Y(ATOMK(COUNT))
          YL = Y(ATOML(COUNT))
          YM = Y(ATOMM(COUNT))
          YN = Y(ATOMN(COUNT))
          ZI = Z(ATOMI(COUNT))
          ZJ = Z(ATOMJ(COUNT))
          ZK = Z(ATOMK(COUNT))
          ZL = Z(ATOML(COUNT))
          ZM = Z(ATOMM(COUNT))
          ZN = Z(ATOMN(COUNT))
C
          OBSSANI = SANIOBS(COUNT)
          ERRSANI = SANIERR(COUNT)
C
C Calculate phi and its derivatives with respect to x,y,z of atoms:
C i, j, k, l, m, and n. This was done with Mathematica !!!
C
C
        o1=-XI+XJ
        o2=o1**2
        o3=-YI+YJ
        o4=o3**2
        o5=-ZI+ZJ
        o6=o5**2
        o7=o2+o4+o6
        o8=sqrt(o7)
        o9=1/o8
        o10=-XM+XN
        o11=o1*o10
        o12=-YM+YN
        o13=o12*o3
        o14=-ZM+ZN
        o15=o14*o5
        o16=o11+o13+o15
        o17=o10**2
        o18=o12**2
        o19=o14**2
        o20=o17+o18+o19
        o21=sqrt(o20)
        o22=1/o21
        o23=1/o7
        o24=o16**2
        o25=1/o20
        o26=o23*o24*o25
        o27=sqrt(1.d0-o26)
        o28=1/o27
        o29=XM-XN
        o30=o7**(3.d0/2.d0)
        o31=1/o30
        o32=o1*o16*o22*o31
        o33=YM-YN
        o34=o16*o22*o3*o31
        o35=ZM-ZN
        o36=o16*o22*o31*o5
        o37=o20**(3.d0/2.d0)
        o38=1/o37
        o39=o10*o16*o38*o9
        o40=o12*o16*o38*o9
        o41=o14*o16*o38*o9
        o42=-XI+XK
        o43=o42**2
        o44=-YI+YK
        o45=o44**2
        o46=-ZI+ZK
        o47=o46**2
        o48=o43+o45+o47
        o49=sqrt(o48)
        o50=-XI+XL
        o51=o50**2
        o52=-YI+YL
        o53=o52**2
        o54=-ZI+ZL
        o55=o54**2
        o56=o51+o53+o55
        o57=sqrt(o56)
        o58=1/o57
        o59=o10*o42
        o60=o12*o44
        o61=o14*o46
        o62=o59+o60+o61
        o63=1/o62
        o64=o10*o50
        o65=o12*o52
        o66=o14*o54
        o67=o64+o65+o66
        o68=o62**2
        o69=1/o68
        o70=o56**(3.d0/2.d0)
        o71=1/o70
        o72=o49*o50*o63*o67*o71
        o73=1/o49
        o74=o42*o58*o63*o67*o73
        o75=1/o56
        o76=o67**2
        o77=o48*o69*o75*o76
        o78=1/(1.d0+o77)
        o79=o49*o52*o63*o67*o71
        o80=o44*o58*o63*o67*o73
        o81=o49*o54*o63*o67*o71
        o82=o46*o58*o63*o67*o73
C
        THETA = acos(o16*o22*o9)
C
        DTXI = -(o28*(o32+o22*o29*o9))
        DTYI = -(o28*(o34+o22*o33*o9))
        DTZI = -(o28*(o36+o22*o35*o9))
        DTXJ = -(o28*(-o32+o10*o22*o9))
        DTYJ = -(o28*(-o34+o12*o22*o9))
        DTZJ = -(o28*(-o36+o14*o22*o9))
        DTXK = 0
        DTYK = 0
        DTZK = 0
        DTXL = 0
        DTYL = 0
        DTZL = 0
        DTXM = -(o28*(o39+o22*o9*(XI-XJ)))
        DTYM = -(o28*(o40+o22*o9*(YI-YJ)))
        DTZM = -(o28*(o41+o22*o9*(ZI-ZJ)))
        DTXN = -(o28*(-o39+o1*o22*o9))
        DTYN =  -(o28*(-o40+o22*o3*o9))
        DTZN = -(o28*(-o41+o22*o5*o9))
C
        PHI = atan(o49*o58*o63*o67)
C
        DPXI = (o29*o49*o58*o63-o29*o49*o58*o67*o69+o72-o74)*o78
        DPYI = o78*(o33*o49*o58*o63-o33*o49*o58*o67*o69+o79-o80)
        DPZI = o78*(o35*o49*o58*o63-o35*o49*o58*o67*o69+o81-o82)
        DPXJ = 0
        DPYJ = 0
        DPZJ = 0
        DPXK = (-(o10*o49*o58*o67*o69)+o74)*o78
        DPYK = o78*(-(o12*o49*o58*o67*o69)+o80)
        DPZK = o78*(-(o14*o49*o58*o67*o69)+o82)
        DPXL = (o10*o49*o58*o63-o72)*o78
        DPYL = o78*(o12*o49*o58*o63-o79)
        DPZL = o78*(o14*o49*o58*o63-o81)
        DPXM = o78*(-(o49*o58*o67*o69*(XI-XK))+o49*o58*o63*(XI-XL))
        DPYM = o78*(-(o49*o58*o67*o69*(YI-YK))+o49*o58*o63*(YI-YL))
        DPZM = o78*(-(o49*o58*o67*o69*(ZI-ZK))+o49*o58*o63*(ZI-ZL))
        DPXN = (o49*o50*o58*o63-o42*o49*o58*o67*o69)*o78
        DPYN = (o49*o52*o58*o63-o44*o49*o58*o67*o69)*o78
        DPZN = (o49*o54*o58*o63-o46*o49*o58*o67*o69)*o78
C
C
C
C
C calculate energy as a function of delta anis
C
C
C
C Casting the coefficients to relate to the correct definition
C COEF0 = DFS
C COEF2 = B/A   i.e. rhombicity
C
C************************
C IF COEF1 is defined as the minimum observed value
C
C COEF1 = -A - 1.5*A*COEF2 + DFS
C and
C       A = (COEF0 - COEF1)/(1. + (1.5*COEF2))
C************************
C
C************************
C IF COEF1 is defined as the magnitude of the axial component
C
C COEF1 = A
C
C
        A = COEF1
        B = A * COEF2
C
        CALCSANI = COEF0+A*(3.*COS(THETA)*COS(THETA)-1.)+
     &            B*(3./2.)*(SIN(THETA)*SIN(THETA)*COS(2.*PHI))
          IF (WHICH.EQ.'ANALYZE') SANICALC(COUNT) = CALCSANI
          DELTASANI = CALCSANI-OBSSANI
          IF (     (ABS(DELTASANI).LT.ERRSANI)
     &        .AND.(POTENTIAL.EQ.SQUARE)) THEN
              E = 0
              DF1 = 0
              DF2 = 0
          ELSE
           E = K1*(DELTASANI**2)
           DF1 = 2*K1*DELTASANI*(A*(-6.)*SIN(THETA)*COS(THETA)
     &                +B*3.*SIN(THETA)*COS(THETA)*COS(2.*PHI))
           DF2 = 2*K1*DELTASANI*(B*(-3.)*SIN(THETA)*SIN(THETA)*
     &           SIN(2.*PHI))
C
          END IF
C
C
C accumulate energy
C
          ESA = ESA + E
C
C now update forces if in energy & force mode
C
ccc      write (6,'(5i5)') count,atomi(count),atomj(count)
          IF (WHICH.NE.'ANALYZE') THEN
           DX(ATOMI(COUNT))=DX(ATOMI(COUNT))+DF1*DTXI+DF2*DPXI
           DY(ATOMI(COUNT))=DY(ATOMI(COUNT))+DF1*DTYI+DF2*DPYI
           DZ(ATOMI(COUNT))=DZ(ATOMI(COUNT))+DF1*DTZI+DF2*DPZI
           DX(ATOMJ(COUNT))=DX(ATOMJ(COUNT))+DF1*DTXJ+DF2*DPXJ
           DY(ATOMJ(COUNT))=DY(ATOMJ(COUNT))+DF1*DTYJ+DF2*DPYJ
           DZ(ATOMJ(COUNT))=DZ(ATOMJ(COUNT))+DF1*DTZJ+DF2*DPZJ
           DX(ATOMK(COUNT))=DX(ATOMK(COUNT))+DF1*DTXK+DF2*DPXK
           DY(ATOMK(COUNT))=DY(ATOMK(COUNT))+DF1*DTYK+DF2*DPYK
           DZ(ATOMK(COUNT))=DZ(ATOMK(COUNT))+DF1*DTZK+DF2*DPZK
           DX(ATOML(COUNT))=DX(ATOML(COUNT))+DF1*DTXL+DF2*DPXL
           DY(ATOML(COUNT))=DY(ATOML(COUNT))+DF1*DTYL+DF2*DPYL
           DZ(ATOML(COUNT))=DZ(ATOML(COUNT))+DF1*DTZL+DF2*DPZL
           DX(ATOMM(COUNT))=DX(ATOMM(COUNT))+DF1*DTXM+DF2*DPXM
           DY(ATOMM(COUNT))=DY(ATOMM(COUNT))+DF1*DTYM+DF2*DPYM
           DZ(ATOMM(COUNT))=DZ(ATOMM(COUNT))+DF1*DTZM+DF2*DPZM
           DX(ATOMN(COUNT))=DX(ATOMN(COUNT))+DF1*DTXN+DF2*DPXN
           DY(ATOMN(COUNT))=DY(ATOMN(COUNT))+DF1*DTYN+DF2*DPYN
           DZ(ATOMN(COUNT))=DZ(ATOMN(COUNT))+DF1*DTZN+DF2*DPZN
          END IF
      END DO
      RETURN
      END
C
C================
      SUBROUTINE READSANI
C
C reads in coupling constant information
C
C by Nico Tjandra and Marius Clore June 1997
C================
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'susc_anis.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'numbers.inc'
C i/o
C local variables
      INTEGER COUNT, SPTR, OLDCLASS, OLDMAXSANI
      DOUBLE PRECISION K1, CUTOFF, COEF0, COEF1, COEF2
      CHARACTER*4 THENAME
C begin
C
C this is used by READSANI2 to hold the selection
C
      SPTR=ALLHP(INTEG4(NATOM))
C
C reset database only if no couplings have been entered
C ie., this is the first time that SANIsotropy has been called
C
      IF (SANIIPTR.EQ.0) THEN
           CALL SANIDEFAULTS
           CALL ALLOCSANI(0, MAXSANI)
      END IF
C
C now read input
C
      CALL PUSEND('SANI>')
      DO WHILE (.NOT.DONE)
           CALL NEXTWD('SANI>')
           CALL MISCOM('SANI>',USED)
           IF (.NOT.USED) THEN
C
           IF (WD(1:4).EQ.'HELP') THEN
C
              CALL CNSHELP('cns-sanisotropy')
C
C Get class name.  Determine if it's an already-defined class.
C Insert a new class if it's not.
C
           ELSE IF (WD(1:4).EQ.'CLAS') THEN
                OLDCLASS = CURCLASS
                CALL NEXTA4('class name =', THENAME)
                MODE = NEW
                DO COUNT = 1, NCLASSES
                     IF (SANICLASSNAMES(COUNT).EQ.THENAME) THEN
                          MODE = UPDATE
                          CURCLASS = COUNT
                     END IF
                END DO
                IF (MODE.EQ.NEW) THEN
C
C make sure you can't add more than the maximum
C number of classes
C
                     IF (OLDCLASS.EQ.MAXSANICLASSES) THEN
                          CALL DSPERR('SANI','Too many classes.')
                          CALL DSPERR('SANI',
     &                      'Increase MAXSANICLASSES and recompile.')
                          CALL WRNDIE(-5, 'READSANI',
     &                      'Too many Anisotropy classes.')
                     END IF
                     NCLASSES = NCLASSES + 1
                     CURCLASS = NCLASSES
                     SANICLASSNAMES(CURCLASS) = THENAME
C
C If this isn't the first class, close off the old class
C
                     IF (NCLASSES.GT.1) THEN
                            SANIASSNDX(OLDCLASS) = NSANI
                     END IF
               END IF
C
C
C set force constant for current class
C
           ELSE IF (WD(1:4).EQ.'FORC') THEN
                CALL NEXTF('force constant =', K1)
C
C start a default class if there isn't one defined
C
                IF (CURCLASS.EQ.0) THEN
                     NCLASSES = 1
                     CURCLASS = 1
                END IF
        WRITE(PUNIT, '(A, A, A, F8.3)')
     &            'Setting force const for class ',
     &            SANICLASSNAMES(CURCLASS), ' to ', K1
                SANIFORCES(CURCLASS) = K1
C
C set coefficient constant for current class
C
           ELSE IF (WD(1:4).EQ.'COEF') THEN
                CALL NEXTF('coefficient 0 =', COEF0)
                CALL NEXTF('coefficient 1 =', COEF1)
                CALL NEXTF('coefficient 2 =', COEF2)
C
C start a default class if there isn't one already defined
C
                IF (CURCLASS.EQ.0) THEN
                     NCLASSES = 1
                     CURCLASS = 1
                END IF
                WRITE(PUNIT, '(A,A,A, F8.3, F8.3, F8.3, F8.3)')
     &            'Setting coefficients for class ',
     &            SANICLASSNAMES(CURCLASS), 'to', COEF0, COEF1, COEF2
C
                SANICOEF0(CURCLASS) = COEF0
                SANICOEF1(CURCLASS) = COEF1
                SANICOEF2(CURCLASS) = COEF2
C
C
C reset couplings database
C
           ELSE IF (WD(1:4).EQ.'RESE') THEN
                CALL SANIHP
                CALL SANIINIT
                CALL ALLOCSANI(0, MAXSANI)
C
C set potential type
C
           ELSE IF (WD(1:4).EQ.'POTE') THEN
                CALL NEXTA4('potential type =', THENAME)
                IF (THENAME.EQ.'SQUA') THEN
                     WRITE(PUNIT, '(A)') 'using square well potential.'
                     POTENTIAL = SQUARE
                ELSE IF (THENAME.EQ.'HARM') THEN
                     WRITE(PUNIT, '(A)') 'using harmonic potential.'
                     POTENTIAL = HARMONIC
                ELSE
                  CALL DSPERR('SANI','unknown potential. Using square.')
                  POTENTIAL = SQUARE
                END IF
C
C
C change number of assignment slots
C
           ELSE IF (WD(1:4).EQ.'NRES') THEN
                OLDMAXSANI = MAXSANI
                CALL NEXTI('number of slots =', MAXSANI)
                CALL ALLOCSANI(OLDMAXSANI, MAXSANI)
                WRITE(PUNIT, '(A, I8)')
     &            'SANI: Allocating space for', MAXSANI,
     &            'number of constraints '
C
C
C read in an assignment
C
           ELSE IF (WD(1:4).EQ.'ASSI') THEN
C
C make sure you can't add more coupling assignments
C than you have slots for
C
                IF (NSANI.EQ.MAXSANI) THEN
                     CALL DSPERR('SANI','Too many assignments.')
                     CALL DSPERR('SANI',
     &                    'Increase NREStraints and run again.')
                     CALL WRNDIE(-1,'SANI>',
     &                    'exceeded allocation for SANI restraints')
                END IF
C
C if there isn't a class specified,
C start a default class
C
                IF (CURCLASS.EQ.0) THEN
                     NCLASSES = 1
                     CURCLASS = 1
                END IF
                CALL READSANI2(HEAP(SANIIPTR), HEAP(SANIJPTR),
     &                  HEAP(SANIKPTR), HEAP(SANILPTR),
     &                  HEAP(SANIMPTR),
     &                  HEAP(SANINPTR), HEAP(SANIOBSPTR),
     &                  HEAP(SANIERRPTR), HEAP(SPTR))
C
C print violations
C
           ELSE IF (WD(1:4).EQ.'PRIN') THEN
                CALL NEXTWD('PRINt>')
                IF (WD(1:4).NE.'THRE') THEN
                  CALL DSPERR('SANI',
     &                        'print expects THREshold parameter.')
                ELSE
                     CALL NEXTF('THREshold =', CUTOFF)
                     IF (CUTOFF.LT.ZERO) THEN
                       CALL DSPERR('SANI', 'cutoff must be positive.')
                       CUTOFF = ABS(CUTOFF)
                     END IF
                     CALL NEXTA4('ALL or CLASs>', THENAME)
                     IF (THENAME(1:3).EQ.'ALL') THEN
                          DO COUNT = 1,NCLASSES
                               PRINTCLASS(COUNT) = .TRUE.
                          END DO
                     ELSE IF (THENAME(1:4).EQ.'CLAS') THEN
                          CALL NEXTA4('class name =', THENAME)
                          DO COUNT = 1,NCLASSES
                               IF (SANICLASSNAMES(COUNT).EQ.
     &                              THENAME) THEN
                                    PRINTCLASS(COUNT) = .TRUE.
                               ELSE
                                    PRINTCLASS(COUNT) = .FALSE.
                               END IF
                          END DO
                     ELSE
                          DO COUNT = 1,NCLASSES
                               PRINTCLASS(COUNT) = .TRUE.
                          END DO
                     END IF
C
           CALL PRINTSANI(CUTOFF, HEAP(CALCSANIPTR), HEAP(SANIOBSPTR),
     &                 HEAP(SANIIPTR), HEAP(SANIJPTR), HEAP(SANIKPTR),
     &                 HEAP(SANILPTR), HEAP(SANIMPTR), HEAP(SANINPTR))
                END IF
C
C check for END statement
C
           ELSE
                CALL CHKEND('SANI>', DONE)
           END IF
           END IF
      END DO
      DONE = .FALSE.
      CALL FREHP(SPTR,INTEG4(NATOM))
      RETURN
      END
C===============
      SUBROUTINE ALLOCSANI (OLDSIZE, NEWSIZE)
C
C resets coupling constant arrays to hold SIZE entries
C
C by Nico Tjandra and Marius Clore June 1997
C===============
      IMPLICIT NONE
C include files
      INCLUDE 'funct.inc'
      INCLUDE 'susc_anis.inc'
C i/o
      INTEGER OLDSIZE, NEWSIZE
C begin
      IF (OLDSIZE.NE.0) THEN
           CALL FREHP(SANIIPTR, INTEG4(OLDSIZE))
           CALL FREHP(SANIJPTR, INTEG4(OLDSIZE))
           CALL FREHP(SANIKPTR, INTEG4(OLDSIZE))
           CALL FREHP(SANILPTR, INTEG4(OLDSIZE))
           CALL FREHP(SANIMPTR, INTEG4(OLDSIZE))
           CALL FREHP(SANINPTR, INTEG4(OLDSIZE))
           CALL FREHP(SANIOBSPTR, IREAL8(OLDSIZE))
           CALL FREHP(SANIERRPTR, IREAL8(OLDSIZE))
           CALL FREHP(CALCSANIPTR, IREAL8(OLDSIZE))
      END IF
C
      SANIIPTR = 0
      SANIJPTR = 0
      SANIKPTR = 0
      SANILPTR = 0
      SANIMPTR = 0
      SANINPTR = 0
      SANIOBSPTR = 0
      SANIERRPTR = 0
      CALCSANIPTR = 0
C
      IF (NEWSIZE.NE.0) THEN
        SANIIPTR = ALLHP(INTEG4(NEWSIZE))
        SANIJPTR = ALLHP(INTEG4(NEWSIZE))
        SANIKPTR = ALLHP(INTEG4(NEWSIZE))
        SANILPTR = ALLHP(INTEG4(NEWSIZE))
        SANIMPTR = ALLHP(INTEG4(NEWSIZE))
        SANINPTR = ALLHP(INTEG4(NEWSIZE))
        SANIOBSPTR = ALLHP(IREAL8(NEWSIZE))
        SANIERRPTR = ALLHP(IREAL8(NEWSIZE))
        CALCSANIPTR = ALLHP(IREAL8(NEWSIZE))
      END IF
C
      RETURN
      END
C==============
      SUBROUTINE SANIDEFAULTS
C
C sets up defaults
C
C by Nico Tjandra and Marius Clore June 1997
C==============
      IMPLICIT NONE
C include files
      INCLUDE 'susc_anis.inc'
C local variables
      INTEGER COUNT
C begin
      SANIIPTR = 0
      MODE = NEW
      MAXSANI = 200
      NSANI = 0
      NCLASSES = 0
      CURCLASS = 0
      POTENTIAL = SQUARE
      DO COUNT = 1, MAXSANICLASSES
           SANICLASSNAMES(COUNT) = 'DEFAULT'
           SANIASSNDX(COUNT) = 0
           SANIFORCES (COUNT) = 50
           SANICOEF0(COUNT) = 1.0
           SANICOEF1(COUNT) = 1.0
           SANICOEF2(COUNT) = 1.0
      END DO
      RETURN
      END
C==============
      SUBROUTINE READSANI2 (ATOMI, ATOMJ, ATOMK, ATOML,
     &         ATOMM, ATOMN, SANIOBS, SANIERR, SEL)
C
C reads actual J coupling assignments into arrays
C
C by Nico Tjandra and Marius Clore June 1997
C==============
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'susc_anis.inc'
      INCLUDE 'mtf.inc'
C i/o
      INTEGER ATOMI(*), ATOMJ(*), ATOMK(*), ATOML(*)
      INTEGER ATOMM(*), ATOMN(*),  SEL(*)
      DOUBLE PRECISION SANIOBS(*), SANIERR(*)
C local variables
      INTEGER NSEL, INSERTPOS, COUNT, CURSTOP, OTHERSTOP, NFLAG
      DOUBLE PRECISION SANIO, SANIE
C begin
C
C if we're in update mode, make a space for the new line
C
      NFLAG = 0
C
      IF (MODE.EQ.UPDATE) THEN
           DO COUNT = NSANI+1, SANIASSNDX(CURCLASS)+1, -1
                ATOMI(COUNT) = ATOMI(COUNT-1)
                ATOMJ(COUNT) = ATOMJ(COUNT-1)
                ATOMK(COUNT) = ATOMK(COUNT-1)
                ATOML(COUNT) = ATOML(COUNT-1)
                ATOMM(COUNT) = ATOMM(COUNT-1)
                ATOMN(COUNT) = ATOMN(COUNT-1)
                SANIOBS(COUNT) = SANIOBS(COUNT-1)
                SANIERR(COUNT) = SANIERR(COUNT-1)
           END DO
           CURSTOP = SANIASSNDX(CURCLASS)
           DO COUNT = 1, NCLASSES
                OTHERSTOP = SANIASSNDX(COUNT)
                IF (OTHERSTOP.GT.CURSTOP) THEN
                     SANIASSNDX(COUNT) = OTHERSTOP + 1
                END IF
           END DO
           SANIASSNDX(CURCLASS) = CURSTOP + 1
           INSERTPOS = CURSTOP
           NSANI = NSANI + 1
      ELSE
           NSANI = NSANI + 1
           INSERTPOS = NSANI
           SANIASSNDX(CURCLASS) = INSERTPOS
      END IF
C
      CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.)
      IF (NSEL.GT.1) THEN
       CALL DSPERR('SANI',
     &             'more than 1 atom in sel for atom i. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
       CALL DSPERR('SANI','Error with atom selection')
       NFLAG = 1
      ENDIF
      CALL MAKIND(SEL, NATOM, NSEL)
      ATOMI(INSERTPOS) = SEL(1)
C
      CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.)
      IF (NSEL.GT.1) THEN
       CALL DSPERR('SANI',
     &             'more than 1 atom in sel for atom j. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
       CALL DSPERR('SANI','Error with atom selection')
       NFLAG = 1
      ENDIF
      CALL MAKIND(SEL, NATOM, NSEL)
      ATOMJ(INSERTPOS) = SEL(1)
C
      CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.)
      IF (NSEL.GT.1) THEN
       CALL DSPERR('SANI',
     &             'more than 1 atom in sel for atom k. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
       CALL DSPERR('SANI','Error with atom selection')
       NFLAG = 1
      ENDIF
      CALL MAKIND(SEL, NATOM, NSEL)
      ATOMK(INSERTPOS) = SEL(1)
C
      CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.)
      IF (NSEL.GT.1) THEN
       CALL DSPERR('SANI',
     &             'more than 1 atom in sel for atom l. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
       CALL DSPERR('SANI','Error with atom selection')
       NFLAG = 1
      ENDIF
      CALL MAKIND(SEL, NATOM, NSEL)
      ATOML(INSERTPOS) = SEL(1)
C
      CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.)
      IF (NSEL.GT.1) THEN
       CALL DSPERR('SANI',
     &             'more than 1 atom in sel for atom l. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
       CALL DSPERR('SANI','Error with atom selection')
       NFLAG = 1
      ENDIF
      CALL MAKIND(SEL, NATOM, NSEL)
      ATOMM(INSERTPOS) = SEL(1)
C
      CALL SELCTA(SEL, NSEL, X, Y, Z,.TRUE.)
      IF (NSEL.GT.1) THEN
       CALL DSPERR('SANI',
     &             'more than 1 atom in sel for atom l. Using first')
      END IF
      IF (NSEL.EQ.0) THEN
       CALL DSPERR('SANI','Error with atom selection')
       NFLAG = 1
      ENDIF
      CALL MAKIND(SEL, NATOM, NSEL)
      ATOMN(INSERTPOS) = SEL(1)
C
      CALL NEXTF('observed Anis =', SANIO)
      CALL NEXTF('error in Anis =', SANIE)
      SANIOBS(INSERTPOS) =SANIO
      SANIERR(INSERTPOS) = SANIE
C
C Check for error atom selection. If there is one than reset the
C counter for restraint
C
      IF (NFLAG.EQ.1) THEN
       NSANI = NSANI-1
      ENDIF
C
      RETURN
      END
C===============
      SUBROUTINE SANIINIT
C
C initializes J-coupling stuff
C
C by nico Tjandra and Marius Clore June 1997
C================
      IMPLICIT NONE
C include files
      INCLUDE 'susc_anis.inc'
C begin
      CALL SANIDEFAULTS
C
      RETURN
      END
C===============
      SUBROUTINE SANIHP
C
C deallocates susceptibility anisotropy stuff
C
C by Gregory Warren 1/9/98
C================
      IMPLICIT NONE
C include files
      INCLUDE 'susc_anis.inc'
C begin
      IF (SANIIPTR.NE.0) THEN
      CALL ALLOCSANI(MAXSANI,0)
      END IF
      RETURN
      END
C=================
      SUBROUTINE PRINTSANI (CUTOFF, SANICALC, SANIOBS,
     &       ATOMI, ATOMJ, ATOMK, ATOML, ATOMM, ATOMN)
C
C prints anisotropy with deltaJ greater than cutoff
C calculates RMS deviation and puts it into $RESULT
C
C by Nico Tjandra and Marius Clore June 1997
C=================
      IMPLICIT NONE
C include files
      INCLUDE 'cns.inc'
      INCLUDE 'susc_anis.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'numbers.inc'
C i/o
      DOUBLE PRECISION CUTOFF, SANICALC(*), SANIOBS(*)
      INTEGER ATOMI(*), ATOMJ(*), ATOMK(*), ATOML(*)
      INTEGER ATOMM(*), ATOMN(*)
C local variables
      DOUBLE PRECISION CALCSANI, OBSSANI, DELTASANI, SANIENERGY
      INTEGER COUNT, CLASS, I, J, K, L, M, N, NUM
      DOUBLE PRECISION RMS, VIOLS, DBPREC
      DOUBLE COMPLEX DUMMY2
      LOGICAL PRINTTHISCLASS
C begin
      RMS = ZERO
      VIOLS = ZERO
      NUM = 0
C
C make sure that the calcJ array is up to date
C
      CALL ESANI(SANIENERGY, 'ANALYZE')
      WRITE (PUNIT, '(A)') 'The following anisotropies have'
      WRITE (PUNIT, '(A)') 'delta SANI greater than the cutoff:'
      WRITE (PUNIT, '(A)') '(calc Anis) (obs Anis) (delta Anis)'
C
C write out first class heading
C
      CLASS = 1
      PRINTTHISCLASS = PRINTCLASS(CLASS)
      IF (PRINTTHISCLASS) THEN
           WRITE (PUNIT, '(A, A)') 'class ', SANICLASSNAMES(1)
      END IF
C
C for every anisotropy entry,
C
      COUNT = 0
      DO WHILE(COUNT.LT.NSANI)
          COUNT=COUNT+1
C
C is this the start of a new class?
C
          IF (SANIASSNDX(CLASS).LT.COUNT) THEN
               CLASS = CLASS + 1
               PRINTTHISCLASS = PRINTCLASS(CLASS)
               IF (SANIASSNDX(CLASS).EQ.SANIASSNDX(CLASS-1)) THEN
                 COUNT = COUNT - 1
               ENDIF
               IF (PRINTTHISCLASS) THEN
                 WRITE (PUNIT, '(A, A)') 'class ', SANICLASSNAMES(CLASS)
               END IF
          END IF
C
C If this assignment is in a class to be printed
C and make sure there is an entry for that class
C
         IF((PRINTTHISCLASS).AND.(SANIASSNDX(CLASS).NE.
     &    SANIASSNDX(CLASS-1))) THEN
C
C always update RMS delta J
C
          CALCSANI = SANICALC(COUNT)
          OBSSANI = SANIOBS(COUNT)
          DELTASANI = (CALCSANI-OBSSANI)
          RMS = RMS + DELTASANI**2
          NUM = NUM + 1
C
C print out delta Js greater than cutoff
C and update number of violations
C
          IF (ABS(DELTASANI).GT.CUTOFF) THEN
               I = ATOMI(COUNT)
               J = ATOMJ(COUNT)
               K = ATOMK(COUNT)
               L = ATOML(COUNT)
               M = ATOMM(COUNT)
               N = ATOMN(COUNT)
               WRITE(PUNIT,'(5X,8(1X,A), 3(F8.3))')
     &               SEGID(M),RESID(M),RES(M),TYPE(M),
     &               SEGID(N),RESID(N),RES(N),TYPE(N),
     &               CALCSANI, OBSSANI, DELTASANI
               VIOLS = VIOLS + ONE
          END IF
       END IF
      END DO

      IF (NUM.GT.0) THEN
          RMS = SQRT(RMS / NUM)
      ELSE
          RMS=0.0
      ENDIF
      DBPREC = NUM
      CALL DECLAR('NUMBER', 'DP', ' ', DUMMY2, DBPREC)
      CALL DECLAR('RESULT', 'DP', ' ', DUMMY2, RMS)
      CALL DECLAR('VIOLATIONS', 'DP', ' ', DUMMY2, VIOLS)
      RETURN
      END
C
      SUBROUTINE SCRASANI
C
C Author: Gregory Warren
      IMPLICIT NONE
C I/O
      INCLUDE 'susc_anis.inc'
C reset diffusion anisotropy database
      IF (SANIIPTR.NE.0) THEN
      WRITE(6,'(A)')
     & ' SCRATC-warning: susceptibility anistropy database erased.'
      CALL SANIHP
      CALL SANIINIT
      END IF
      RETURN
      END
C
