$DEBUG C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C NSYM - PREDICTIVE PROGRAM FOR SYMMETRIC TOPS WITH ONE QUDRUPOLAR NUCLEUS C AND SPIN-ROTATION COUPLING C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C THIS PROGRAM WILL PREDICT THE SYMMETRIC ROTOR SPECTRUM C UP TO J+1--J 10--9 C (FOR EXTENDED VERSION SEE JAMAN.FTN) C C This version is practically identical to the original NSYM of the Exeter C Microwave Group, the only changes are some tidying up and adaptation C for running on a PC as well as on a mainframe C C C Version 3.1/89 ZBIGNIEW KISIEL C C C NOTE: Tests carried out while producing this documentation showed that C NSYM will work satisfactorily for I=3/2 (CH3Cl,CH3Br) but WILL NOT C satisfactorily reproduce results for I=5/2 (and possibly other C values of I). Sensible reproduction of results could not be C obtained for: C 1. 13.CD3I - Mallinson JMS 55,94(1975) (without sextic cd and sr) C 2. CH3I - Burie et al Mol.Phys. 32,289(1976) (with sextic cd and sr) C C In addition the sign of spin-rotation terms in the program has C been reversed to conform to general usage. The FORTRAN coding C also leaves some scope for improvement... C C Z.K. C____________________________________________________________________________ C C RUNNING ON A MAINFRAME: C 1. Remove $DEBUG if present - if the MSF version is being converted C 2. Insert the acceptable edit descriptor for stopping the cursor C in all FORMAT statements using it. In MS-Fortran this is '\', C on DEC and many other machines it is ',$'. This is not C standardised in FORTRAN-77. C RUNNING ON A PC (MICROSOFT FORTRAN): C 1. Replace all ',$' (or other) in FORMAT statements with '\' C C_____________________________________________________________________________ C C THE INPUT SHOULD BE FREE FORMAT AS FOLLOWS: C C line 1 TITLE C line 2 JHD : set to 1 if HJ,HJK&HKJ required, otherwise to 0 C line 3 JSR : set to 1 if CN&CK required, otherwise to 0 C line 4 UMX : maximum transition freq. for calculation C line 4 B DJ DJK (HJ HJK HKJ CN CK) X.o I KMAX C C - Terms in brackets in line 4 are to be included as determined by C values of JHD and JSR C - X.o is eQq C - Nuclear spin I should be in units of 1/2 C - KMAX is the maximum value of K for the calculation C - All constants are to be in the same units (MHz) C C The program is semi-interactive in that it will ask for the input and C output file names and display some intermediate information on the C terminal. C C___________________________________________________________________________ C C IMPLICIT REAL*8 (A-E,G-H,O-Z),INTEGER (F) INTEGER*2 N REAL*8 H(16,16),E(16,26),F,I,J,K REAL*8 T(50),IN(50),ISUM INTEGER FN1(50),FN2(50),FSL CHARACTER*30 FILNAM CHARACTER*72 COM CHARACTER*2 DIVISR DIVISR='/2' C WRITE(*,156) 156 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | N S Y M - Prediction of spectrum of a linear/symmetric', * T79,'|'/ * ' | top with one quadrupolar nucleus', * T79,'|'/ * ' |',76(1H_),'|'/' version 3.1/1989',T64,'Zbigniew KISIEL'/) c 5556 WRITE(*,5555)' INPUT FILE NAME: ' 5555 FORMAT(1X/1X,A\) READ(*,'(A)',ERR=5556)FILNAM OPEN(5,FILE=FILNAM,status='old',ERR=5556) C C...read data C READ(5,'(A)')COM READ(5,*)JHD READ(5,*)JSR READ(5,*)UMX IF(JHD.NE.0 .AND.JSR.NE.0)GO TO 51 IF(JHD.NE.0)GO TO 52 IF(JSR.NE.0)GO TO 53 READ(5,*)B,DJ,DJK,Q,II,KMAX HJ=0.0D0 HJK=0.0D0 HKJ=0.0D0 CN=0.0D0 CK=0.0D0 GO TO 54 C 51 READ(5,*)B,DJ,DJK,HJ,HJK,HKJ,CN,CK,Q,II,KMAX GO TO 54 C 52 READ(5,*)B,DJ,DJK,HJ,HJK,HKJ,Q,II,KMAX CN=0.0D0 CK=0.0D0 GO TO 54 C 53 READ(5,*)B,DJ,DJK,CN,CK,Q,II,KMAX HJ=0.0D0 HJK=0.0D0 HKJ=0.0D0 C C 54 WRITE(*,5555)'OUTPUT FILE NAME: ' READ(*,'(A)',ERR=54)FILNAM OPEN(6,FILE=FILNAM,STATUS='UNKNOWN',ERR=54) c WRITE(6,190)COM 190 FORMAT(1X,72(1H_)//1X,A/1X,72(1H_)//' Prediction of symmetric', 1 ' top nuclear quadrupole pattern for:'//) JMAX=10 JJMAX=JMAX*2 FMIN=0 IF(((-1)**II).LT.0)FMIN=1 IF(FMIN.EQ.0)DIVISR=' ' ISPIN=II IF(FMIN.EQ.0)ISPIN=II/2 WRITE(6,7000)ISPIN,DIVISR,B,DJ,DJK,HJ,HJK,HKJ,CN,CK,Q 7000 FORMAT(4X,'I =',I2,A2//4X,'B =',F14.6,' MHz' 1 /4X,'DJ =',3PF14.6,' kHz'/4X,'DJK =',F14.6,' kHz' 2 /4X,'HJ =',6PF14.6,' Hz'/4X,'HJK =',F14.6,' Hz' 3 /4X,'HKJ =',F14.6,' Hz'/4X,'CN =',3PF14.6,' kHz' 4 /4X,'CK =',F14.6,' kHz'/4X,'eQq =',0PF14.6,' MHz') C C write(6,6000) write(*,'(//'' Now working on K = ''\)') DO 60 KK=0,KMAX write(*,'(i5\)')kk DO 2000 L=1,16 DO 2000 M=1,24 E(L,M)=0.0D0 2000 CONTINUE DO 10 FF=FMIN,JJMAX+II,2 DO 1000 L=1,16 DO 1000 M=1,16 H(L,M)=0.0D0 1000 CONTINUE JL=IABS(FF-II)/2 JH=(FF+II)/2 IF(KK.GT.JL)JL=KK IF(KK.GT.JH)GOTO 10 DO 20 JJ=JL,JH F=DBLE(FLOAT(FF))/2.0D0 I=DBLE(FLOAT(II))/2.0D0 J=JJ K=KK JK=JJ+1 C=F*(F+1.0D0)-I*(I+1.0D0)-J*(J+1.0D0) IF(J.NE.0.0D0) GO TO 11 H(JK,JK)=-Q*(0.75D0*C*(C+1.0D0)-I*(I+1.0D0)* 1 J*(J+1.0D0))/(2.0D0*I*(2.0D0*I-1.0D0)* 1 (2.0D0*J-1.0D0)*(2.0D0*J+3.0D0)) GO TO 12 11 H(JK,JK)=Q*(3.0D0*K**2/(J*(J+1.0D0))-1.0D0) 1 *(0.75D0*C*(C+1.0D0)-I*(I+1.0D0)*J*(J+1.0D0)) 1 /(2.0D0*I*(2.0D0*I-1.0D0)*(2.0D0*J-1.0D0) 1 *(2.0D0*J+3.0D0)) 12 CONTINUE IF(JK+1.GT.12)GOTO 90 IF(J.NE.0.0D0) GO TO 13 H(JK,JK)=0.0D0 GOTO 14 13 H(JK,JK+1)=(3.0D0*Q*K*(F*(F+1.0D0)-I*(I+1.0D0) 1 -J*(J+2.0D0)) 1 /(8.0D0*I*(2.0D0*I-1.0D0)*J*(J+2.0D0))) 1 *DSQRT((1.0D0-K**2/(J+1.0D0)**2)*(I+J+F+2.0D0) 1 *(J+F-I+1.0D0)*(I+F-J)*(I+J-F+1.0D0) 1 /((2.0D0*J+1.0D0)*(2.0D0*J+3.0D0))) 14 CONTINUE 90 IF(JK+2.GT.12)GOTO 100 H(JK,JK+2)=(3.0D0*Q/(16.0D0*I*(2.0D0*I-1.0D0) 1 *(2.0D0*J+3.0D0))) 1 *DSQRT((1.0D0-K**2/(J+1.0D0)**2) 2 *(1.0D0-K**2/(J+2.0D0)**2)*(I+J+F+3.0D0) 1 *(I+J+F+2.0D0)*(I+J-F+2.0D0)*(I+J-F+1.0D0) 3 *(J+F-I+2.0D0)*(J+F-I+1.0D0)*(I+F-J) 4 *(I+F-J-1.0D0)/((2.D0*J+1.D0)*(2.D0*J+5.D0))) 100 H(JK,JK)=H(JK,JK)+B*J*(J+1.0D0)-DJ*J**2 1 *(J+1.0D0)**2-DJK*J*(J+1.0D0)*K**2 IF(JHD.NE.0 .AND. JSR.NE.0) GO TO 71 IF(JHD.NE.0) GO TO 72 IF(JSR.NE.0) GO TO 73 GO TO 20 71 IF(J.EQ.0.0D0) GO TO 74 H(JK,JK)=H(JK,JK)+HJ*J**3*(J+1.0D0)**3 1 +HJK*J**2*(J+1.0D0)**2*K**2 1 +HKJ*J*(J+1.0D0)*K**4-0.5D0* 2 (CN+((CK-CN)*K*K)/(J*(J+1.0D0)))*C GO TO 20 74 H(JK,JK)=H(JK,JK)+HJ*J**3*(J+1.0D0)**3 1 +HJK*J**2*(J+1.0D0)**2*K**2 1 +HKJ*J*(J+1.0D0)*K**4-0.5D0*CN*C GO TO 20 72 H(JK,JK)=H(JK,JK)+HJ*J**3*(J+1.0D0)**3 1 +HJK*J**2*(J+1.0D0)**2*K**2 1 +HKJ*J*(J+1.0D0)*K**4 GO TO 20 73 IF(J.EQ.0.0D0) GO TO 75 H(JK,JK)=H(JK,JK)-0.5D0*(CN+(CK-CN*K*K)/(J*(J+1.0D0)))*C GO TO 20 75 H(JK,JK)=H(JK,JK)-0.5D0*CN*C 20 CONTINUE DO 70 MM=1,16 IF(MM+1.GT.12)GOTO 80 H(MM+1,MM)=H(MM,MM+1) 80 IF(MM+2.GT.12)GOTO 70 H(MM+2,MM)=H(MM,MM+2) 70 CONTINUE N=16 CALL ADIAGJ(H,N) FK=FF+1 DO 30 JK=JL+1,JH+1 E(JK,FK)=H(JK,JK) C WRITE(6,*)JK,FK,E(JK,FK) 30 CONTINUE 10 CONTINUE ASSIGN 3000 TO FSL MF=1 IF(FMIN.EQ.1)ASSIGN 4000 TO FSL IF(FMIN.EQ.1)MF=2 DO 40 JJ=KK,JMAX-1 J2=JJ*2 FL1=IABS(J2-II) FU1=J2+II J1=(JJ+1)*2 FL2=IABS(J1-II) FU2=J1+II ISUM=0.0D0 N=0 DO 110 F1=FL1,FU1,2 DO 110 F2=FL2,FU2,2 IF(IABS(F2-F1)/2.GT.1)GOTO 110 N=N+1 T(N)=E(JJ+2,F2+1)-E(JJ+1,F1+1) J=JJ J=J+1.0D0 F=F1/2.0D0 IF((F2-F1)/2)120,130,140 120 F=F-1.0D0 IN(N)=(J-F+I)*(J-F+I-1.0D0)*(J-F-I-1.0D0) 1 *(J-F-I-2.0D0)/(F+1.0D0) GOTO 150 130 IN(N)=-1.0D0*(J+F+I+1.0D0)*(J+F-I)*(J-F+I)*(J-F-I-1.0D0) 1 *(2.0D0*F+1.0D0)/(F*(F+1.0D0)) GOTO 150 140 F=F+1.0D0 IN(N)=(J+F+I+1.0D0)*(J+F+I)*(J+F-I)*(J+F-I-1.0D0)/F 150 ISUM=ISUM+IN(N) 110 CONTINUE N=0 DO 50 F1=FL1,FU1,2 DO 50 F2=FL2,FU2,2 IF(IABS(F2-F1)/2.GT.1)GOTO 50 N=N+1 IN(N)=(IN(N)/ISUM) FF1=F1 FF2=F2 IF(FMIN.EQ.0)FF1=FF1/2 IF(FMIN.EQ.0)FF2=FF2/2 FN1(N)=FF1 FN2(N)=FF2 50 CONTINUE NSUM=N 170 IFLAG=0 DO 160 N=1,NSUM-1 IF(T(N+1).GE.T(N))GOTO 160 IFLAG=1 TEMP=T(N) T(N)=T(N+1) T(N+1)=TEMP TEMP=IN(N) IN(N)=IN(N+1) IN(N+1)=TEMP TEMP=FN1(N) FN1(N)=FN1(N+1) FN1(N+1)=TEMP TEMP=FN2(N) FN2(N)=FN2(N+1) FN2(N+1)=TEMP 160 CONTINUE IF(IFLAG.EQ.1)GOTO 170 IF(T(1).GT.UMX) GO TO 40 WRITE(6,'(1x)') DO 41 NN=1,NSUM JJP1=JJ+1 INTERM=(FN2(NN)-FN1(NN))/MF WRITE(6,FSL)JJP1,KK,FN2(NN),JJ,KK,FN1(NN),T(NN),IN(NN), 1 INTERM 41 CONTINUE 40 CONTINUE 60 CONTINUE C 3000 FORMAT(3X,I2,2X,I2,2X,I2,6X,I2,2X,I2,2X,I2,7X,F13.5,7X,F7.5, 1 8X,I2) 4000 FORMAT(3X,I2,2X,I2,2X,I2,'/2',4X,I2,2X,I2,2X,I2,'/2',5X,F13.5, 1 7X,F7.5,8X,I2) 6000 FORMAT(//4X,'J''',2X,'K''',2X,'F''',6X,'J',3X,'K',3X,'F',7X, 1 'FREQUENCY/MHz',5X,'INTENSITY',5X,6H(F'-F),/) C CLOSE(5) write(6,'(1x,72(1h_))') CLOSE(6) WRITE(*,'(1x///)') C STOP END c c____________________________________________________________________________ c SUBROUTINE INV(AA,N,D,IAA,JAA) REAL*8 AA(IAA,JAA),A(100) INTEGER*2 L(10),M(10) REAL*8 D,BIGA,HOLD,DABS C IJ=0 DO 5000 J=1,N DO 5000 I=1,N IJ=IJ+1 A(IJ)=AA(I,J) 5000 CONTINUE C C SEARCH FOR LARGEST ELEMENT C D=1.0 NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 20 J=K,N IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I 10 IF(DABS(BIGA)-DABS(A(IJ)))15,20,20 15 BIGA=A(IJ) L(K)=I M(K)=J 20 CONTINUE C C INTERCHANGE ROWS C J=L(K) IF(J-K)35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) A(JI)=HOLD 30 CONTINUE C C INTERCHANGE COLUMNS C 35 I=M(K) IF(I-K)45,45,38 38 JP=N*(I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) A(JI)=HOLD 40 CONTINUE C C DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS C CONTAINED IN BIGA) C 45 IF(BIGA)48,46,48 46 D=0.0 RETURN 48 DO 55 I=1,N IF(I-K)50,55,50 50 IK=NK+I A(IK)=A(IK)/(-BIGA) 55 CONTINUE C C REDUCE MATRIX C DO 65 I=1,N IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K)60,65,60 60 IF(J-K)62,65,62 62 KJ=IJ-I+K A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE C C DIVIDE ROW BY PIVOT C KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K)70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE C C PRODUCT OF PIVOTS C D=D*BIGA C C REPLACE PIVOT BY RECIPROCAL C A(KK)=1.0/BIGA 80 CONTINUE C C FINAL ROW AND COLUMN INTERCHANGE C K=N 100 K=(K-1) IF(K)150,150,105 105 I=L(K) IF(I-K)120,120,108 108 JQ=N*(K-1) JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) A(JI)=HOLD 110 CONTINUE 120 J=M(K) IF(J-K)100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) A(JI)=HOLD 130 CONTINUE GOTO 100 C 150 DO 5001 I=1,N DO 5001 J=1,N IJ=N*(J-1)+I AA(I,J)=A(IJ) 5001 CONTINUE C RETURN END c c____________________________________________________________________________ c SUBROUTINE ADIAGJ(A,N) C C JACOBI DIAGONALISATION - no eigenvectors, eigenvalues only C IMPLICIT REAL*8(A-H,O-Z),INTEGER*2(I-N) DIMENSION A(N,N) PI=3.141592653589793D0 DO 500 ITRAP=1,500 X=0.0 DO 4 I=1,N DO 301 J=1,N IF(I.EQ.J) GO TO 101 IF(DABS(A(I,J)).GE.DABS(X)) GO TO 100 GO TO 101 100 X=A(I,J) II=I KK=J 101 CONTINUE 301 CONTINUE 4 CONTINUE IF(A(II,II).EQ.A(KK,KK)) GO TO 200 OE=2.0*A(II,KK)/(A(II,II)-A(KK,KK)) PHI=0.50D0*DATAN(OE) GO TO 201 200 IF(A(II,KK).EQ.0.0D0) GO TO 15 IF(A(II,KK).GT.0.0) PHI=PI/4.0D0 IF(A(II,KK).LT.0.0) PHI=-1.0*PI/4.0D0 201 R=(A(II,II)-A(KK,KK))**2+4.0D0*A(II,KK)**2 R=DSQRT(R) IF(A(II,II).GE.A(KK,KK)) SIGMA=1.0D0 IF(A(II,II).LT.A(KK,KK)) SIGMA=-1.0D0 D1=0.50D0*(A(II,II)+A(KK,KK)+SIGMA*R) D2=0.50D0*(A(II,II)+A(KK,KK)-SIGMA*R) SINE=DSIN(PHI) COSI=DCOS(PHI) DO 700 I=1,N A1=A(II,I) A2=A(KK,I) A(II,I)=A1*COSI+A2*SINE A(I,II)=A(II,I) A(KK,I)=-A1*SINE+A2*COSI A(I,KK)=A(KK,I) 700 CONTINUE A(II,II)=D1 A(KK,KK)=D2 A(II,KK)=0.0D0 A(KK,II)=0.0D0 IF(DABS(X).LE.1.0D-10) GO TO 15 500 CONTINUE WRITE(6,501) 501 FORMAT('Algorithim fails to converge ') 15 RETURN END c____________________________________________________________________________ c____________________________________________________________________________