PROGRAM INS
C     LAST CHANGES 3-19-96 PETR BOUR                                  
      IMPLICIT REAL*8 (A-H,O-Z)                                       
      PARAMETER (N0=1300,NIT=20000)
      REAL*8 SV(N0,3*N0,3),R(3,N0),FRQ(3*N0),Z(N0),
     1SG(N0),AI(3,3),TI(3),AM(N0),SGC(N0),AINS(NIT),OT(NIT)
      INTEGER A,B                                                     
      CHARACTER*4 ATID
C
      BOHR=0.52917715                                                 
      Z0=0.0
C
      OPEN(3,FILE='INS.OUT')
      WRITE(6,6000)
      WRITE(3,6000)
6000  FORMAT(
     1' Approximate intensities of inelastic neutron scattering',/,/,
     2'          Input:    F.INP',/,
     3'                 FTRY.INP',/,/,
     4'         Output:  INS.TAB',/,
     5'                  INS.OUT',/)
C
      OPEN(17,FILE='F.INP',STATUS='OLD')                       
      READ(17,*)NQ,N3,NOAT
      DO 35 K=1,NOAT                                                  
      READ(17,*)IDEN,(R(I,K),I=1,3)                                   
      DO 34 I=1,3                                                     
34    R(I,K)=R(I,K)/BOHR                                              
      Z(K)=Z0
      IF (IDEN.EQ.6)Z(K)=6.0D0                                        
      IF (IDEN.EQ.8)Z(K)=8.0D0                                        
      IF (IDEN.EQ.7)Z(K)=7.0D0                                        
      IF (IDEN.EQ.1)Z(K)=1.0D0                                        
      IF (IDEN.EQ.16)Z(K)=16.0D0                                       
      IF (IDEN.EQ.15)Z(K)=15.0D0                                       
      IF (IDEN.EQ.9)Z(K)=9.0D0                                       
      IF (IDEN.EQ.17)Z(K)=17.0D0
      IF (DABS(Z(K)).LT.1.0D-10)THEN    
       WRITE(6,*)' UNKNOWN ATOM ',IDEN,' PROGRAM STOPPED'             
       STOP                                                           
      ENDIF                                                           
35    CONTINUE                                                        
      WRITE(*,*)' Geometry read in, ',NOAT,' atoms'
      WRITE(3,*)' Geometry read in, ',NOAT,' atoms'
      READ(17,*)                                                      
      DO 45 LL=1,NOAT                                                 
      DO 45 I=1,NQ                                                    
45    READ (17,*) IAT,IMD,SV(IAT,I,1),SV(IAT,I,2),SV(IAT,I,3)         
      READ(17,*)                                                      
      READ(17,1717)(FRQ(I),I=1,NQ)                                    
1717  FORMAT(6F11.3)                                                  
      WRITE(3,*)' S-matrix read in ',NQ,' transitions'
      WRITE(*,*)' S-matrix read in ',NQ,' transitions'
      CLOSE(17)
C
      OPEN(4,FILE='FTRY.INP',STATUS='OLD')
      READ(4,*)
      READ(4,*)NOCS
      DO 1 I=1,NOCS
1     READ(4,*)
      READ(4,4000)(ADUM,I=1,NOAT)
4000  FORMAT(6F12.6)
      READ(4,*)
      READ(4,4000)(AM(I),I=1,NOAT)
      CLOSE(4)
C
      WRITE(3,3000)
3000  FORMAT('     Isotope     INS scattering cross-section',/,
     1       '                 coherent      total    ',/, 
     2       ' --------------------------------------------')
C     "Protein Crystalography", T.L.Blundel, L.N.Johnson, Academic Press,
C     London, 1993; pp468-469.
C     Given in Barns=10^-24 cm^2)
      TOL=0.3
      DO 2 I=1,NOAT
      SG(I)=5.00
      SGC(I)=5.00
      ATID='XXXX'
      IF(ABS(AM(I)-12.0).LT.TOL)THEN
       ATID=' C12'
       SG(I)=5.51
       SGC(I)=5.5
      ENDIF
      IF(ABS(AM(I)-1.0).LT.TOL)THEN
       ATID=' H 1'
       SG(I)=81.5
       SGC(I)=1.79
      ENDIF
      IF(ABS(AM(I)-2.0).LT.TOL)THEN
       ATID=' H 2'
       SG(I)=7.6
       SGC(I)=5.4
      ENDIF
      IF(ABS(AM(I)-13.0).LT.TOL)THEN
       ATID=' C13'
       SG(I)=5.5
       SGC(I)=4.5
      ENDIF
      IF(ABS(AM(I)-14.0).LT.TOL)THEN
       ATID=' N14'
       SG(I)=11.4
       SGC(I)=11.0
      ENDIF
      IF(ABS(AM(I)-16.0).LT.TOL)THEN
       ATID=' O16'
       SG(I)=4.24
       SGC(I)=4.2
C      this was for O15 (?missprint)
      ENDIF
      IF(ABS(AM(I)-31.0).LT.TOL)THEN
       ATID=' P31'
       SG(I)=3.6
       SGC(I)=3.5
      ENDIF
      IF(ABS(AM(I)-32.0).LT.TOL)THEN
       ATID=' S32'
       SG(I)=1.2
       SGC(I)=1.2
      ENDIF
      IF(ABS(AM(I)-35.0).LT.TOL)THEN
       ATID='Cl35'
       SG(I)=15.0
       SGC(I)=12.2
      ENDIF
2     WRITE(3,3001)I,ATID,SGC(I),SG(I)
3001  FORMAT(I4,1X,A4,2F15.2)
      WRITE(3,*)' --------------------------------------------'
      WRITE(6,*)' -- Choose the frequency profile ------------'
      WRITE(6,*)'    1    ...   exponential decay'
      WRITE(6,*)'    2    ...   rational decay   ' 
      WRITE(6,*)' --------------------------------------------'
      READ(*,*)IW
      IF(IW.EQ.1)CALL WE(NQ,FRQ,SV,AINS,NOAT,NST,OT,SG)
      IF(IW.EQ.2)CALL WR(NQ,FRQ,SV,AINS,NOAT,NST,OT,SG)
C
C
      OPEN(4,FILE='INS.TAB')
      WRITE(4,4001)
      WRITE(3,4001)
4001  FORMAT(' Inelastic neutron scattering',/,
     1       '     freq[cm-1]   INS intensity',/,
     2       '-------------------------------')
      AMAX=0.0D0
      DO 4 I=1,NST
4     IF(AMAX.LT.AINS(I))AMAX=AINS(I)
      SF=0.0D0
      IF(AMAX.GT.0.0000000001D0)SF=100.0D0/AMAX
      DO 3 I=1,NST
      WRITE(3,3002)I,OT(I),AINS(I)*SF
3     WRITE(4,3002)I,OT(I),AINS(I)*SF
3002  FORMAT(I4,F12.3,F12.4)
      WRITE(4,4002)
      WRITE(3,4002)
4002  FORMAT('-------------------------------')
      CLOSE(4)
C
      CLOSE(3)
      WRITE(6,*)' >> Good bye <<'
      STOP                                                            
      END                                                             
c     ===================================================================
      SUBROUTINE WE(NQ,FRQ,SV,AINS,NOAT,NST,OT,SG)
      IMPLICIT REAL*8 (A-H,O-Z)                                       
      PARAMETER (N0=1300,NIT=20000)
      DIMENSION AINS(NIT),SV(N0,3*N0,3),FRQ(NQ),OT(NIT),SG(N0)  
c
      WRITE(6,*)' Exponential decay'
      WRITE(6,*)' Input alpha (about 0.06):'
      READ(*,*)ALPHA
      C=1.0D0/16.132D0
c
      NST=NQ
      DO 41 I=1,NQ                                                    
      W=FRQ(I)
      OT(I)=W
      PHI=1.0D0/3.0D0*C*W*EXP(-C*ALPHA*W)
      SK=0.0
      DO 5 L=1,NOAT
      DO 5 iA=1,3
5     SK=SK+SG(L)*SV(L,I,iA)*SV(L,I,iA)
41    AINS(I)=PHI*SK
c
      RETURN
      END
c     ===================================================================
      SUBROUTINE WR(NQ,FRQ,SV,AINS,NOAT,NST,OT,SG)
      IMPLICIT REAL*8 (A-H,O-Z)                                       
      PARAMETER (N0=1300,NEM=5,NIS=20000,NIT=20000)
c
c     NEM ... maximum number of excitation
C     NIS ... maximum number of initial states
C     NIT ... maximum number of transitions
C
      DIMENSION AINS(NIT),SV(N0,3*N0,3),FRQ(NQ),O(NIS),
     1OT(NIT),SS(3*N0),SG(N0)
      INTEGER*1 ST(NIS,3*N0)
      CHARACTER*1 OK
c
      WRITE(6,*)' Rational decay with Boltzman factors'
      WRITE(6,*)' The temperature [K]:'
      READ(*,*)T
c
c     temperature in cm-1:
      TC=T*0.69502D0
c
      REC=4.6D0*TC
      WRITE(6,6000)REC
6000  FORMAT(' Energy cut-off [recommended ',F9.1,' cm-1]: ')
      READ(*,*)WMAX
c
      WRITE(6,*)' Final neutron energy [probably 24 cm-1]:'
      READ(*,*)WFF
c
      WRITE(6,*)' Angle between kf and ki [say 135 deg]:'
      READ(*,*)gg
      cf=cos(gg*3.141592653589D0/180.0)
c
c     Prepare list of states:
      NS=0
      DO 1 I=1,NIS
      O(I)=0.0D0
      DO 1 IM=1,NQ
1     ST(I,IM)=0
c
c     Ground state
      NS=1
c 
c     Monoexcited states
      DO 21 I1=1,NQ
      W1=FRQ(I1)
      IF(W1.LT.WMAX)THEN
       NS=NS+1
       IF(NS.GT.NIS)GOTO 999
       O(NS)=W1
       ST(NS,I1)=1
c
c      Doubly-excited states
       DO 22 I2=I1,NQ
       W2=FRQ(I2)
       W12=W1+W2
       IF(W12.LT.WMAX)THEN
        NS=NS+1
        IF(NS.GT.NIS)GOTO 999
        O(NS)=W12
        ST(NS,I2)=ST(NS,I2)+1
        ST(NS,I1)=ST(NS,I1)+1
c
c       Triple-excited states
        DO 23 I3=I2,NQ
        W3=FRQ(I3)
        W123=W12+W3
        IF(W123.LT.WMAX)THEN
         NS=NS+1
         IF(NS.GT.NIS)GOTO 999
         O(NS)=W123
         ST(NS,I3)=ST(NS,I3)+1
         ST(NS,I2)=ST(NS,I2)+1
         ST(NS,I1)=ST(NS,I1)+1
c
c        Quadruple-excited states
         DO 24 I4=I3,NQ
         W4=FRQ(I4)
         W1234=W123+W4
         IF(W1234.LT.WMAX)THEN
          NS=NS+1
          IF(NS.GT.NIS)GOTO 999
          O(NS)=W1234
          ST(NS,I4)=ST(NS,I4)+1
          ST(NS,I3)=ST(NS,I3)+1
          ST(NS,I2)=ST(NS,I2)+1
          ST(NS,I1)=ST(NS,I1)+1
         ENDIF
24       CONTINUE
        ENDIF
23      CONTINUE
       ENDIF
22     CONTINUE
      ENDIF
21    CONTINUE
c
c     precalculate mode cross sections
      DO 5 I=1,NQ
      SS(I)=0.0D0
      DO 5 L=1,NOAT
      DO 5 IA=1,3
5     SS(I)=SS(I)+SG(L)*SV(L,I,IA)*SV(L,I,IA)
c
      WRITE(3,3002)
3002  FORMAT('transition energy mode Boltzmann',/,80(1H-))
      NST=NQ
      DO 66 I=1,NQ
      OT(I)=FRQ(I)
66    AINS(I)=0.0D0
c
c     Loop over initial states
      DO 3 I=1,NS
      BOLTZ=EXP(-O(I)/TC)
c
c     Loop over allowed transitions,(R)NF - the FINAL quantum number for mode IT
      DO 3 IT=1,NQ
      RNF=REAL(ST(I,IT)+1)
      WF=FRQ(IT)
      AINS(IT)=AINS(IT)+RNF*BOLTZ
c     WRITE(3,3000)NST,WF,IT,BOLTZ
c3     WRITE(3,3001)(ST(I,II),II=1,NQ)
3000  FORMAT(I5,F8.2,I4,F9.6,' from :')
3001  FORMAT(80I1)
3     CONTINUE
c
      DO 77 I=1,NQ
      WF=OT(I)
77    AINS(I)=AINS(I)*(2.0D0*WFF+WF-2.0D0*cf*sqrt(WFF*(WFF+WF)) )
     1           /WF/SQRT(WFF+WF)*SS(I)
c
      RETURN
c
999   WRITE(3,*)' Too many ground states'
      WRITE(6,*)' Too many ground states'
      goto 777
c
888   WRITE(3,*)' Too many transitions'
      WRITE(6,*)' Too many transitions'
777   close(3)
      STOP
c
      END