PROGRAM PARANMRD !dotto con tau-N pesati secondo Boltmann ! versione risc IBM C PARAMAGNETIC ENHANCEMENT IN NMRD PROFILE C C THIS PROGRAM REQUIRES THE INPUT FILE PAR.DAT C 0 BEFORE A PARAMETER MEANS IT HAS TO BE ASSUMED AS CONSTANT C 1 BEFORE A PARAMATER MEANS IT HAS TO BE CHANGED IN THE FITTING C C C INTERNAL FITTING IS PERFORMED IN THE PARAMETERS: C d , Ddiff , RK , A/h , MOLAR FRACTION , TAUM C C SUBROUTINES: C FUNCZFS: SET PARAMETERS IN FITTING PROCEDURE C FUNCINT: SET PARAMETERS IN INTERNAL FITTING PROCEDURE C DIAG: WRITE ENERGY MATRIX AND CALCULATE EXPECTATION VALUES C POWELL: FITTING PROCEDURE C CONNECTED SUBROUTINES: C LINMIN C F1DIM C MNBRAK C BRENT C XISTEP C POWELLINT: INTERNAL FITTING PROCEDURE C F01BCF....X04BAF: CALCULATE EIGENVALUES AND EIGENFUNCTIONS C GAUINT: PERFORME INTEGRATION ON ANGLES C TUNO: PERFORME CALCULATION OF T1M C TDUE: PERFORME CALCULATION OF T2M C TUNOISO: PERFORME CALCULATION OF T1M IN AXIAL CASE C C THETA AND PHI ARE THE POLAR ANGLES DEFINING THE WATER PROTON DIRECTION C IN THE MOLECULAR FRAME C C BETA AND GAMMA ARE THE EULER ANGLES DEFINING THE MOLECULAR FRAME WITH C RESPECT TO THE LAB FRAME IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10) COMMON /SET/SET COMMON /RK10/ SPIN, SI COMMON /GAMMAH/ GAMMAI COMMON /B31/ TPUNO(500) COMMON /B32/ PP(500),Z(500) COMMON /TAUDELTA/ TAUDELTA COMMON /B4/ NVMEM,NPT(10),NPTOT COMMON /WATER/ ACQ COMMON /T1T2/ IREL COMMON /INDEX/ INDEX COMMON /ALFASTEP/ ALFASTEP COMMON /STAMPA/INDEXSTAMPA COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND COMMON/INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10) & ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10) COMMON /C1M/DM(10),DDM(10),CONCM(10) COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10), & TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10), & GXM(10),GYM(10),GZM(10) COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10), & TAUMM(10,10) COMMON /MOLFRAZM/ AMOLFRAM(10) COMMON /CONTATM/ ACONTM(10) COMMON /CICLE/ NVEST real*4 b,DIF,xMSAT,TAUN,HANI,TEMP1 common /super/ b,DIF,xMSAT,TAUN,HANI,TEMP1 common /super2/ omega0,ppar,qpar real*4 omega0,ppar,qpar COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10), & B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10), & B16(10),B17(10),B18(10),B19(10),B20(10),B21(10) CHARACTER*20 FILENAME COMMON/TEMPERATURE/ TEMP(10) COMMON /TMSTART/ TM11(10),TM21(10) C DIMENSION=MAX NUMBER OF PARAMETERS (21) DIMENSION P(21) DIMENSION P1(21) COMMON /PPAR/ P2(21) DIMENSION XI(21,21) INDEX=1 INDEXSTAMPA=0 C CONSTANTS READ IN FILE PAR.DAT OPEN (1, STATUS = 'OLD', FILE = 'PARC.DAT') OPEN (4, FILE="PAR.OUT") C OUTPUT FILE READ(1,'(A)')FILENAME C NUCLEAR MOLECULAR SPIN READ(1,*)SI C GAMMA OF THE INVESTIGATED PARTICLE READ(1,*)GAMMAI C ELECTRON SPIN READ(1,*)SPIN C T1 OR T2 CALCULATION READ(1,*)IREL C LIMITS OF THE FIELD READ(1,*)X1,X2,X3 IF(X3.EQ.1)THEN XMIN=X1 XMAX=X2 ELSE XMIN=LOG10(GAMMAI*X1/6.283) XMAX=LOG10(GAMMAI*X2/6.283) ENDIF C NUMBER OF POINTS TO BE CALCULATED READ(1,*)NUMPUN IF(XMIN.EQ.XMAX)NUMPUN=1 C NUMBER OF SETS OF DATA FOR FITTING READ(1,*)SET IF(SET.EQ.0) SET=1 C TEMPERATURE READ(1,*)(TEMP(K),K=1,SET) J=1 IND=1 IND1=1 IND2=1 NV=0 NVEST=0 NPLUS=0 NPLUS2=0 C CORRELATION TIMES READ(1,*)B1(J), (TAUS0M(J,K),K=1,2),TAUDELTA IF(B1(J).GE.2)THEN TS1=TAUS0M(J,1) TS2=TAUS0M(J,2) DO K=1,SET TAUS0M(J,K)=TS1*EXP(TS2/TEMP(K)) END DO IF(B1(J).EQ.3)NPLUS=NPLUS+1 ENDIF IF(B1(J).EQ.1)THEN P(IND)=TAUS0M(J,1) P1(IND1)=TAUS0M(J,1) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF IF(B1(J).EQ.2)THEN P(IND)=TS1 P1(IND1)=TS1 P(IND+1)=TS2 P1(IND1+1)=TS2 WRITE(4,'(2X,4(E10.4,2X))') TS1,TS2 IND=IND+2 IND1=IND1+2 ENDIF READ(1,*)B2(J), (TAURM(J,K),K=1,2) IF(B2(J).GE.2)THEN TR1=TAURM(J,1) TR2=TAURM(J,2) DO K=1,SET TAURM(J,K)=TR1*EXP(TR2/TEMP(K)) !Stokes END DO IF(B2(J).EQ.3) NPLUS=NPLUS+1 ENDIF IF(B2(J).EQ.1)THEN P(IND)=TAURM(J,1) P1(IND1)=TAURM(J,1) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF IF(B2(J).EQ.2)THEN P(IND)=TR1 P1(IND1)=TR1 P(IND+1)=TR2 P1(IND1+1)=TR2 WRITE(4,'(2X,4(E10.4,2X))') TR1,TR2 IND=IND+2 IND1=IND1+2 ENDIF READ(1,*)B3(J), (TAUVM(J,K),K=1,2) IF(B3(J).GE.2)THEN TV1=TAUVM(J,1) TV2=TAUVM(J,2) DO K=1,SET TAUVM(J,K)=TV1*EXP(TV2/TEMP(K)) END DO IF(B3(J).EQ.3) NPLUS=NPLUS+1 ENDIF IF(B3(J).EQ.1)THEN P(IND)=TAUVM(J,1) P1(IND1)=TAUVM(J,1) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF IF(B3(J).EQ.2)THEN P(IND)=TV1 P1(IND1)=TV1 P(IND+1)=TV2 P1(IND1+1)=TV2 WRITE(4,'(2X,4(E10.4,2X))') TV1,TV2 IND=IND+2 IND1=IND1+2 ENDIF IF (TAURM(J,1).EQ.0.AND.TAUVM(J,1).EQ.0)THEN TAUCM(J)=TAUS0M(J,1) IF(B1(J).EQ.1) P(IND-1)=TAUS0M(J,K) IF(B1(J).EQ.1) P1(IND1-1)=TAUS0M(J,K) ENDIF C PARAMETERS OF ZERO FIELD SPLITTING READ(1,*)B11(J), DPARAM(J) DPARAM(J) = PI2*VL*DPARAM(J) IF(B11(J).EQ.1)THEN P(IND)=DPARAM(J) P1(IND1)=DPARAM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2 INDDPARA(J)=IND IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B12(J), EPARAM(J) EPARAM(J) = PI2*VL*EPARAM(J) IF(B12(J).EQ.1)THEN P(IND)=EPARAM(J) P1(IND1)=EPARAM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2 INDEPARA(J)=IND IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B19(J), S4M(J) S4M(J) = PI2*VL*S4M(J) IF(B19(J).EQ.1)THEN P(IND)=S4M(J) P1(IND1)=S4M(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2 INDS4(J)=IND IND=IND+1 IND1=IND1+1 ENDIF C PARAMETER OF G-TENSOR READ(1,*)B17(J), GSER GXM(J)=GSER/2.003 IF(B17(J).EQ.1)THEN P(IND)=GXM(J) P1(IND1)=GXM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B18(J), GSER GYM(J)=GSER/2.003 IF(B18(J).EQ.1)THEN P(IND)=GYM(J) P1(IND1)=GYM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B20(J), GSER GZM(J)=GSER/2.003 IF(B20(J).EQ.1)THEN P(IND)=GZM(J) P1(IND1)=GZM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF C PARAMETERS OF HYPERFINE COUPLING READ(1,*)B9(J), AXM(J) C CONVERTION FROM CM-1 TO S-1.RAD AXM(J)=PI2*VL*AXM(J) IF(B9(J).EQ.1)THEN P(IND)=AXM(J) P1(IND1)=AXM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2 INDAPAR(J)=IND IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B10(J), AYM(J) AYM(J)=PI2*VL*AYM(J) IF(B10(J).EQ.1)THEN P(IND)=AYM(J) P1(IND1)=AYM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2 INDAPER(J)=IND IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B21(J), AZM(J) AZM(J)=PI2*VL*AZM(J) IF(B21(J).EQ.1)THEN P(IND)=AZM(J) P1(IND1)=AZM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/VL/PI2 INDAPER2(J)=IND IND=IND+1 IND1=IND1+1 ENDIF C PARAMETERS OF OUTER-SPHERE READ(1,*)B13(J), DM(J),super,dm2,dm3 DM(J)=DM(J)*1.E-8 IF(B13(J).EQ.1)THEN P(IND)=DM(J) P2(IND2)=DM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)/1.E-8 INDED(J)=IND IND=IND+1 IND2=IND2+1 ENDIF if(super.eq.3)then IF(B13(J).EQ.2)THEN P(IND)=dm2 IND=IND+1 ENDIF IF(B13(J).EQ.3)THEN P(IND)=dm3 IND=IND+1 ENDIF IF(B13(J).EQ.4)THEN P(IND)=dm(1) IND=IND+1 P(IND)=dm2 IND=IND+1 P(IND)=dm3 IND=IND+1 ENDIF IF(B13(J).EQ.5)THEN P(IND)=dm2 IND=IND+1 P(IND)=dm3 IND=IND+1 ENDIF endif READ(1,*)B14(J), DDM(J),dd2,dd3 IF(B14(J).EQ.1)THEN P(IND)=DDM(J) P2(IND2)=DDM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND2=IND2+1 ENDIF if(super.ge.2)then IF(B14(J).EQ.2)THEN P(IND)=dd2 IND=IND+1 ENDIF IF(B14(J).EQ.3)THEN P(IND)=dd3 IND=IND+1 ENDIF IF(B14(J).EQ.4)THEN P(IND)=DDM(1) IND=IND+1 P(IND)=dd2 IND=IND+1 P(IND)=dd3 IND=IND+1 ENDIF IF(B14(J).EQ.5)THEN P(IND)=dd2 IND=IND+1 P(IND)=dd3 IND=IND+1 ENDIF IF(B14(J).EQ.6)THEN P(IND)=DDM(1) IND=IND+1 P(IND)=dd3 IND=IND+1 ENDIF IF(B14(J).EQ.7)THEN P(IND)=DDM(1) IND=IND+1 P(IND)=dd2 IND=IND+1 ENDIF endif READ(1,*)B15(J), CONCM(J) IF(B15(J).EQ.1)THEN P(IND)=CONCM(J) P2(IND2)=CONCM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND2=IND2+1 ENDIF C NUMBER OF DIFFERENT SITES READ(1,*) ACQ C PARAMETERS FOR DIFFERENT SITES DO J=1,ACQ READ(1,*)B4(J), (TAUMM(J,K),K=1,2) IF(B4(J).GE.2)THEN TM1=TAUMM(J,1) TM2=TAUMM(J,2) TM11(J)=TAUMM(J,1) TM21(J)=TAUMM(J,2) DO K=1,SET TAUMM(J,K)=TM1*EXP(TM2/TEMP(K)) END DO IF(B4(J).EQ.3) NPLUS2=1 ENDIF IF(B4(J).EQ.1)THEN P(IND)=TAUMM(J,1) P2(IND2)=TAUMM(J,1) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND2=IND2+1 ENDIF IF(B4(J).EQ.2)THEN P(IND)=TM1 P2(IND2)=TM1 P(IND+1)=TM2 P2(IND2+1)=TM2 WRITE(4,'(2X,4(E10.4,2X))') TM1,TM2 IND=IND+2 IND2=IND2+2 ENDIF READ(1,*)B5(J), AMOLFRAM(J) AMOLFRAM(J)=AMOLFRAM(J)*1.E-3/111. IF(B5(J).EQ.1)THEN P(IND)=AMOLFRAM(J) P2(IND2)=AMOLFRAM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND)*111./1.E-3 INDAMOLFRA(J)=IND IND=IND+1 IND2=IND2+1 ENDIF READ(1,*)B6(J), RKM(J) IF(B6(J).EQ.1)THEN P(IND)=RKM(J) P2(IND2)=RKM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND2=IND2+1 ENDIF READ(1,*)B16(J), ACONTM(J) IF(B16(J).EQ.1)THEN P(IND)=ACONTM(J) P2(IND2)=ACONTM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND2=IND2+1 ENDIF READ(1,*)B7(J), THETAM(J) IF(B7(J).EQ.1)THEN P(IND)=THETAM(J) P1(IND1)=THETAM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF READ(1,*)B8(J), PHIM(J) IF(B8(J).EQ.1)THEN P(IND)=PHIM(J) P1(IND1)=PHIM(J) WRITE(4,'(2X,4(E10.4,2X))') P(IND) IND=IND+1 IND1=IND1+1 ENDIF C NUMBER OF FITTING PARAMETERS NV=NV+B1(J)+B2(J)+B3(J)+B4(J)+B5(J)+B6(J)+B7(J)+B8(J)+B9(J)+ & B10(J)+B11(J)+B12(J)+B13(J)+B14(J)+B15(J)+B16(J)+B17(J) & +B18(J)+B19(J)+B20(J)+B21(J)-NPLUS*3-NPLUS2*3 C NUMBER OF FITTING PARAMETERS IN EXTERNAL CICLE NVEST=NVEST+B1(J)+B2(J)+B3(J)+B7(J)+B8(J)+B9(J)+ & B10(J)+B11(J)+B12(J)+B17(J) & +B18(J)+B19(J)+B20(J)+B21(J)-NPLUS*3 NPLUS=0 NPLUS2=0 END DO NVMEM=NV DPARATOT=0. EPARATOT=0. APERTOT=0. APERTOT2=0. APARTOT=0. ACONTOT=0. ACONIND=0. DO J=1,ACQ DPARATOT=DPARAM(J)+DPARATOT EPARATOT=EPARAM(J)+EPARATOT APERTOT=AZM(J)+APERTOT APERTOT2=AXM(J)+APERTOT2 APARTOT=AYM(J)+APARTOT ACONTOT=ACONTM(J)+ACONTOT END DO C DEFINITION OF NMX: DIMENSION OF ENERGY MATRIX IF(DPARATOT.EQ.0..AND.EPARATOT.EQ.0. & AND.GX.EQ.GZ.AND.GX.EQ.GY.AND.SPIN.EQ.0.5)THEN NMX=2.*(2*SI+1.) ELSE IF(APERTOT.EQ.0.AND.APARTOT.EQ.0.AND.APERTOT2.EQ.0.AND. & GX.EQ.GZ.AND.GX.EQ.GY.AND.EPARATOT.EQ.0)THEN SI=0.5 NMX=2*SPIN+1. ELSE NMX = (2*SI + 1)*(2*SPIN+1) ENDIF ENDIF IF(IREL.NE.1) NMX = (2*SI + 1)*(2*SPIN+1) IF(ACONTOT.NE.0) THEN IF(APERTOT.NE.0.OR.APARTOT.NE.0.OR.APERTOT2.NE.0.OR. & GX.NE.GZ.OR.GX.NE.GY.OR.EPARATOT.NE.0.OR.DPARATOT.NE.0)THEN NMX = (2*SI + 1)*(2*SPIN+1) ACONIND=1. ENDIF ENDIF READ(1,*) (NPT(K),K=1, SET) READ(1,*) FTOL READ(1,*) ALFASTEP C Z: FREQUENCIES, PP: RATE NPTOT=0 DO K=1,SET NPTOT=NPTOT+NPT(K) END DO DO 11 I=1,NPTOT READ(1,*) Z(I),PP(I) 11 CONTINUE CLOSE(1) OO=10 b=dm(1) DIF=ddm(1) xMSAT=dd2 TAUN=taus0m(1,1) HANI=dd3 TEMP1=temp(1) ppar=dm2 qpar=dm3 omega0=hani IF (NPT(1).EQ.0)GOTO 250 C STARTING FITTING PROCEDURE if(super.eq.2.and.npt(1).ne.0)then if(b14(1).eq.1.or.b14(1).eq.2.or.b14(1).EQ.3)B14C=1 if(b14(1).eq.5.or.b14(1).eq.6.or.b14(1).EQ.7)B14C=2 if(b14(1).eq.4)B14C=3 if(b14(1).eq.0)B14C=0 np=b1(1)+b13(1)+b14c N=np ITER=1000 CALL suPOWELL(P,XI,N,NP,FTOL,ITER,FRET,NMX) write(6,*)'Parameters:' do i=1,np write(6,*)p(i) end do goto 1234 endif if(super.eq.3.and.npt(1).ne.0)then if(b13(1).eq.1.or.b13(1).eq.2.or.b13(1).EQ.3)B13B=1 if(b13(1).eq.4)B13B=3 if(b13(1).eq.5)B13B=2 if(b14(1).eq.1.or.b14(1).eq.2.or.b14(1).EQ.3)B14C=1 if(b14(1).eq.5.or.b14(1).eq.6.or.b14(1).EQ.7)B14C=2 if(b14(1).eq.4)B14C=3 if(b14(1).eq.0)B14C=0 np=b1(1)+b14c+b13b N=np ITER=1000 CALL suPOWELL(P,XI,N,NP,FTOL,ITER,FRET,NMX) write(6,*)'Parameters:' do i=1,np write(6,*)p(i) end do goto 1234 endif IF(NVEST.EQ.0)THEN CALL FUNCZFS(P2,FUNC,NMX,NV) NP=NV N=NV ITER=1000 CALL POWELLINT(P2,XI,N,NP,FTOL,ITER,FRET,NMX) DO J=1,NP WRITE(6,'(2X,4(E10.4,2X))') P2(J) END DO ELSE NP=NVEST N=NVEST ITER=1000 CALL POWELL(P1,XI,N,NP,FTOL,ITER,FRET,NMX) ENDIF WRITE(4,*) 'ERROR=', FRET/(NPTOT-NV) DO J=1,ACQ IF(B5(J).EQ.1)P2(INDAMOLFRA(J))=P2(INDAMOLFRA(J))*111/1.E-3 IF(B9(J).EQ.1)P1(INDAPAR(J))=P1(INDAPAR(J))/PI2/VL IF(B10(J).EQ.1)P1(INDAPER(J))=P1(INDAPER(J))/PI2/VL IF(B21(J).EQ.1)P1(INDAPER2(J))=P1(INDAPER2(J))/PI2/VL IF(B11(J).EQ.1)P1(INDDPARA(J))=P1(INDDPARA(J))/PI2/VL IF(B12(J).EQ.1)P1(INDEPARA(J))=P1(INDEPARA(J))/PI2/VL IF(B13(J).EQ.1)P2(INDED(J))=P2(INDED(J))/1.E-8 IF(B19(J).EQ.1)P1(INDS4(J))=P1(INDS4(J))/PI2/VL END DO C WRITE RESULTS OF FITTING PROCEDURE IF(NVEST.NE.0)THEN JI=1 IF(B1(1).EQ.2)THEN DO K=1,SET WRITE(4,'(2X,4(E10.4,2X))')P1(JI)*EXP(P1(JI+1)/TEMP(K)) END DO JI=JI+2 ENDIF IF(B2(1).EQ.2)THEN DO K=1,SET WRITE(4,'(2X,4(E10.4,2X))')P1(JI)*EXP(P1(JI+1)/TEMP(K)) END DO JI=JI+2 ENDIF IF(B3(1).EQ.2)THEN DO K=1,SET WRITE(4,'(2X,4(E10.4,2X))')P1(JI)*EXP(P1(JI+1)/TEMP(K)) END DO JI=JI+2 ENDIF DO J=JI,NVEST WRITE(4,'(2X,4(E10.4,2X))') P1(J) END DO ENDIF DO JJ=1,ACQ IF(NV-NVEST.NE.0)THEN JI=1 IF(B4(JJ).EQ.2)THEN DO K=1,SET WRITE(4,'(2X,4(E10.4,2X))')P2(JI)*EXP(P2(JI+1)/TEMP(K)) END DO JI=JI+2 ENDIF DO J=JI,NV-NVEST WRITE(4,'(2X,4(E10.4,2X))') P2(J) END DO ENDIF END DO WRITE(4,*) 'MAGN. FIELD, OBSED., CAL. ' DO 1 I=1,NPTOT WRITE(4,'(2X,3(F8.3,2X))') Z(I),PP(I)/CONCM(1)*0.001, & TPUNO(I)/CONCM(1)*0.001 1 CONTINUE CLOSE(4) DO J=1,ACQ IF(B5(J).EQ.1)P2(INDAMOLFRA(J))=P2(INDAMOLFRA(J))/111*1.E-3 IF(B9(J).EQ.1)P1(INDAPAR(J))=P1(INDAPAR(J))*PI2*VL IF(B10(J).EQ.1)P1(INDAPER(J))=P1(INDAPER(J))*PI2*VL IF(B21(J).EQ.1)P1(INDAPER2(J))=P1(INDAPER2(J))*PI2*VL IF(B11(J).EQ.1)P1(INDDPARA(J))=P1(INDDPARA(J))*PI2*VL IF(B12(J).EQ.1)P1(INDEPARA(J))=P1(INDEPARA(J))*PI2*VL IF(B13(J).EQ.1)P2(INDED(J))=P2(INDED(J))*1.E-8 IF(B19(J).EQ.1)P1(INDS4(J))=P1(INDS4(J))*PI2*VL END DO OO=0 1234 continue 250 CONTINUE C CALCULATION OF THE CURVE OPEN (44, FILE=FILENAME) DO K=1,SET NPT(K)=NUMPUN END DO INDEXSTAMPA=1 DO K=1,SET DO I=1,NPT(K) ZE = XMIN + (XMAX - XMIN)*FLOAT(I)/FLOAT(NPT(K)) ADD=0 DO IJK=1,K-1 ADD=ADD+NPT(IJK) END DO PP(I+1+ADD)=PP(I+ADD) Z(I+ADD) = 10.**ZE/1000000. END DO END DO NVEST=30+OO if(super.eq.1)CALL FUNCZFS(P1,FUNC,NMX,NV) if(super.ge.2)write(6,*)'Superparamagnetic routine' if(super.eq.3)CALL FUNCsuper(P,FUNC,NMX,NV) if(super.eq.2)CALL FUNCsuper2(P,FUNC,NMX,NV) CLOSE(44) STOP END SUBROUTINE FUNCZFS(P,FUNC,NMX,NV) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(NV) DIMENSION XI2(21,21) PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10) COMMON /SET/SET COMMON /PPAR/ P2(21) COMMON /RK10/ SPIN, SI COMMON /WATER/ ACQ COMMON /STAMPA/INDEXSTAMPA COMMON /STEPGAMMA/ STEPGAMMA COMMON /B4/ NVMEM,NPT(10),NPTOT COMMON /B31/ TPUNO(500) COMMON /B32/ PP(500),Z(500) COMMON /C1M/DM(10),DDM(10),CONCM(10) COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10), & TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),GXM(10), & GYM(10),GZM(10) COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10), & TAUMM(10,10) COMMON /MOLFRAZM/ AMOLFRAM(10) COMMON /CONTATM/ ACONTM(10) COMMON /INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10) & ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10) COMMON /C1/D,DD,CONC COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4 COMMON /GTENSOR/ GX,GY,GZ COMMON /TAU1/ TAUS0 COMMON /TAU/ TAUR,TAUV,TAUM COMMON /MOLFRAZ/ AMOLFRA COMMON /CONTAT/ ACONT COMMON /GAMMAH/ GAMMAI COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS COMMON /CICLE/ NVEST COMMON /TMAT/ TMAT(500,10),TMATCONT(500,10),TMATCROSS(500,10) COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10), & B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10), & B16(10),B17(10),B18(10),B19(10),B20(10),B21(10) COMMON/TEMPERATURE/ TEMP(10) DIMENSION TS1(10),TS2(10),TR1(10),TR2(10),TV1(10),TV2(10) COMMON /TMSTART/ TM11(10),TM21(10) C SET PARAMETERS FB=0. FBW=0. IF(NVEST.LT.30)THEN WRITE(6,*)'PARAMETERS: ' DO J=1,ACQ IF(B5(J).EQ.1)P(INDAMOLFRA(J))=P(INDAMOLFRA(J))*111/1.E-3 IF(B9(J).EQ.1)P(INDAPAR(J))=P(INDAPAR(J))/PI2/VL IF(B10(J).EQ.1)P(INDAPER(J))=P(INDAPER(J))/PI2/VL IF(B21(J).EQ.1)P(INDAPER2(J))=P(INDAPER2(J))/PI2/VL IF(B11(J).EQ.1)P(INDDPARA(J))=P(INDDPARA(J))/PI2/VL IF(B12(J).EQ.1)P(INDEPARA(J))=P(INDEPARA(J))/PI2/VL IF(B13(J).EQ.1)P(INDED(J))=P(INDED(J))/1.E-8 IF(B19(J).EQ.1)P(INDS4(J))=P(INDS4(J))/PI2/VL END DO DO I=1,NV WRITE(6,'(2X,E10.4)')P(I) END DO DO J=1,ACQ IF(B5(J).EQ.1)P(INDAMOLFRA(J))=P(INDAMOLFRA(J))/111*1.E-3 IF(B9(J).EQ.1)P(INDAPAR(J))=P(INDAPAR(J))*PI2*VL IF(B10(J).EQ.1)P(INDAPER(J))=P(INDAPER(J))*PI2*VL IF(B21(J).EQ.1)P(INDAPER2(J))=P(INDAPER2(J))*PI2*VL IF(B11(J).EQ.1)P(INDDPARA(J))=P(INDDPARA(J))*PI2*VL IF(B12(J).EQ.1)P(INDEPARA(J))=P(INDEPARA(J))*PI2*VL IF(B13(J).EQ.1)P(INDED(J))=P(INDED(J))*1.E-8 IF(B19(J).EQ.1)P(INDS4(J))=P(INDS4(J))*PI2*VL END DO ENDIF DO 223 K=1,SET DO I=1, NPT(K) IND=1 TPUNOTOT=0 J=1 TAUS0=TAUS0M(J,1) TAUR=TAURM(J,1) TAUV=TAUVM(J,1) IF(TAUR.EQ.0.AND.TAUV.EQ.0.)TAUC=TAUCM(J) AX=AXM(J) AY=AYM(J) AZ=AZM(J) DPARA=DPARAM(J) EPARA=EPARAM(J) GX=GXM(J) GY=GYM(J) GZ=GZM(J) S4=S4M(J) IF(INDEXSTAMPA.EQ.0.OR.NVEST.EQ.40)THEN D=DM(J) DD=DDM(J) CONC=CONCM(J) ENDIF C PARAMETERS************************************************************** IF(B1(J).EQ.1)THEN TAUS0=P(IND) IND=IND+1 ENDIF IF(B1(J).EQ.2)THEN TS1(K)=P(IND) TS2(K)=P(IND+1) IND=IND+2 ENDIF IF(B2(J).EQ.1)THEN TAUR=P(IND) IND=IND+1 ENDIF IF(B2(J).EQ.2)THEN TR1(K)=P(IND) TR2(K)=P(IND+1) IND=IND+2 ENDIF IF(B3(J).EQ.1)THEN TAUV=P(IND) IND=IND+1 ENDIF IF(B3(J).EQ.2)THEN TV1(K)=P(IND) TV2(K)=P(IND+1) IND=IND+2 ENDIF IF(B1(J).EQ.1.AND.TAUR.EQ.0.AND.TAUV.EQ.0)THEN TAUC=P(IND-1) ENDIF IF(B11(J).EQ.1)THEN DPARA=P(IND) IND=IND+1 ENDIF IF(B12(J).EQ.1)THEN EPARA=P(IND) IND=IND+1 ENDIF IF(B19(J).EQ.1)THEN S4=P(IND) IND=IND+1 ENDIF IF(B17(J).EQ.1)THEN GX=P(IND) IND=IND+1 ENDIF IF(B18(J).EQ.1)THEN GY=P(IND) IND=IND+1 ENDIF IF(B20(J).EQ.1)THEN GZ=P(IND) IND=IND+1 ENDIF IF(B9(J).EQ.1)THEN AX=P(IND) IND=IND+1 ENDIF IF(B10(J).EQ.1)THEN AY=P(IND) IND=IND+1 ENDIF IF(B21(J).EQ.1)THEN AZ=P(IND) IND=IND+1 ENDIF IND2=1 !!!!!!!!!!!!!!!!!!!!!!!!! IF(B13(J).EQ.1) IND2=IND2+1 IF(B14(J).EQ.1) IND2=IND2+1 IF(B15(J).EQ.1) IND2=IND2+1 !!!!!!!!!!!!!!!!! DO J=1,ACQ IF(INDEXSTAMPA.EQ.0.OR.NVEST.EQ.40)THEN TAUM=TAUMM(J,1) AMOLFRA=AMOLFRAM(J) RK=RKM(J) ACONT=ACONTM(J) ENDIF IF(INDEXSTAMPA.EQ.1)THEN !!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(B4(J).EQ.0)TAUM=TAUMM(J,1) IF(B4(J).EQ.1)THEN TAUM=P2(IND2) IND2=IND2+1 ENDIF IF(B4(J).EQ.2)THEN TM11(J)=P2(IND2) TM21(J)=P2(IND2+1) IND2=IND2+2 ENDIF IF(B5(J).EQ.0)AMOLFRA=AMOLFRAM(J) IF(B5(J).EQ.1)THEN AMOLFRA=P2(IND2) IND2=IND2+1 ENDIF IF(B6(J).EQ.0)RK=RKM(J) IF(B6(J).EQ.1)THEN RK=P2(IND2) IND2=IND2+1 ENDIF IF(B16(J).EQ.0)ACONT=ACONTM(J) IF(B16(J).EQ.1)THEN ACONT=P2(IND2) IND2=IND2+1 ENDIF ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!! THETA=THETAM(J) PHI=PHIM(J) IF(B7(J).EQ.1)THEN THETA=P(IND) IND=IND+1 ENDIF IF(B8(J).EQ.1)THEN PHI=P(IND) IND=IND+1 ENDIF IF(B1(J).EQ.2)TAUS0=TS1(1)*EXP(TS2(1)/TEMP(K)) IF(B2(J).EQ.2)TAUR=TR1(1)*EXP(TR2(1)/TEMP(K)) !stokes IF(B3(J).EQ.2)TAUV=TV1(1)*EXP(TV2(1)/TEMP(K)) IF(B4(J).EQ.2)TAUM=TM11(J)*EXP(TM21(J)/TEMP(K)) IF(B1(J).EQ.3)TAUS0=TAUS0M(J,K) IF(B2(J).EQ.3)TAUR=TAURM(J,K) IF(B3(J).EQ.3)TAUV=TAUVM(J,K) IF(B4(J).EQ.3)TAUM=TAUMM(J,K) C******************************************************************************* ADD=0 DO IJK=1,K-1 ADD=ADD+NPT(IJK) END DO IPLUS=I+ADD C PROTON LARMOR FREQUENCY BZ=Z(IPLUS)*1000000 C CONSTANTS IN DIPOLAR RELAXATION RK1=10.*2.46502D-52/(2.6752D8)**2*GAMMAI**2*(RK*1.E-10)**(-6) IF (AMOLFRA.EQ.0.)RK1=10./(SPIN*(SPIN+1.)*2./15.)/TAUS0*1.E-9 IF(RK.EQ.0) RK1=0. C CONSTANT OF CONTACT RELAXATION CONA=(ACONT*6.28*1.E6*1.0546E-34) IF(TAUC.EQ.0.AND.TAUS0.EQ.0) GOTO 56 CALL GAUINT (BZ,TAUM,NMX) C STORE CONTRIBUTIONS FOR TUNO TMAT(IPLUS,J)=TMUNO TMATCONT(IPLUS,J)=TMUNOCONT TMATCROSS(IPLUS,J)=TMUNOCROSS C CALCULATION OF TPUNO TMUNO=TMUNO*RK1+CONA*CONA*TMUNOCONT+SQRT(RK1)*CONA*TMUNOCROSS TPUNO1=1./(1./TMUNO+TAUM)*AMOLFRA IF(AMOLFRA.EQ.0.)TPUNO1=1./(1./TMUNO+TAUM) TPUNO(IPLUS)=TPUNO1/CONC*0.001 TPUNOTOT=TPUNOTOT+TPUNO(IPLUS) END DO 56 CONTINUE C CALCULATION OF OUTER-SPHERE CONTRIBUTION TERM=0 IF(DD.NE.0)THEN V=BZ*2.*3.14 TAUD=D**2/DD CZ=SQRT(2*V*TAUD) ZZ=SQRT(2*V*TAUD*658.) GEI=(1.+5.*CZ/8.+CZ**2/8.)/(1.+CZ+CZ**2/2.+CZ**3/6.+4.* & CZ**4/81.+CZ**5/81.+CZ**6/648.) GES=(1.+5.*ZZ/8.+ZZ**2/8.)/(1.+ZZ+ZZ**2/2.+ZZ**3/6.+4.* & ZZ**4/81.+ZZ**5/81.+ZZ**6/648.) PRIMO=32.*3.14/405.*(2.6752E4*1.7608E7*1.0546E-27)**2*6.022E20 SECONDO=SPIN*(SPIN+1)*CONC/(D*DD) TERZO=(3.*GEI+7.*GES) TERM=(PRIMO*SECONDO*TERZO) ENDIF TPUNOTOT=TPUNOTOT+TERM TPUNO(IPLUS)=TPUNOTOT C DIFFERENCE BETWEEN EXPERIMENTAL AND FITTING VALUES FB=((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FB FBW=SQRT((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FBW IF (STEPGAMMA.NE.1) WRITE(6,'(2X,2(E10.4))')Z(IPLUS),TPUNO(IPLUS) IF(INDEXSTAMPA.EQ.1)WRITE(44,'(2X,2(E10.4))')Z(IPLUS),TPUNO(IPLUS) END DO 223 CONTINUE FUNC=FBW IF(INDEXSTAMPA.EQ.0)WRITE(6,*)' ERROR', & '=',FBW/(NPTOT-NVMEM),'**********' IF(NV.NE.NVMEM)THEN C INTERNAL FITTING PROCEDURE NP2=NVMEM-NV N2=NVMEM-NV ITER2=1000 FRET2=0. CALL POWELLINT(P2,XI2,N2,NP2,FTOL,ITER2,FRET2,NMX) WRITE(6,*)' ERROR', & '=', FRET2/(NPTOT-NVMEM) DO J=1,NP2 WRITE(6,'(2X,E10.4)')P2(J) END DO FUNC=FRET2 ENDIF RETURN END SUBROUTINE FUNCINT(P,FUNC,NMX,NV) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(NV) PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10) COMMON /SET/SET COMMON /RK10/ SPIN, SI COMMON /WATER/ ACQ COMMON /STAMPA/INDEXSTAMPA COMMON /STEPGAMMA/ STEPGAMMA COMMON /B4/ NVMEM,NPT(10),NPTOT COMMON /B31/ TPUNO(500) COMMON /B32/ PP(500),Z(500) COMMON /C1M/DM(10),DDM(10),CONCM(10) COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10), & TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),GXM(10), & GYM(10),GZM(10) COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10), & TAUMM(10,10) COMMON /MOLFRAZM/ AMOLFRAM(10) COMMON /CONTATM/ ACONTM(10) COMMON/INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10) & ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10) COMMON /C1/D,DD,CONC COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4 COMMON /GTENSOR/ GX,GY,GZ COMMON /TAU1/ TAUS0 COMMON /TAU/ TAUR,TAUV,TAUM COMMON /MOLFRAZ/ AMOLFRA COMMON /CONTAT/ ACONT COMMON /GAMMAH/ GAMMAI COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS COMMON /CICLE/ NVEST COMMON /TMAT/ TMAT(500,10),TMATCONT(500,10),TMATCROSS(500,10) COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10), & B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10), & B16(10),B17(10),B18(10),B19(10),B20(10),B21(10) COMMON/TEMPERATURE/ TEMP(10) DIMENSION TM1(10),TM2(10) C SET PARAMETERS FB=0. FBW=0. DO 223 K=1,SET DO I=1, NPT(K) IND=1 TPUNOTOT=0 J=1 C PARAMETERS OF INTERNAL FITTING********************************************* IF(B13(J).EQ.1)THEN D=P(IND) IND=IND+1 ELSE D=DM(J) ENDIF IF(B14(J).EQ.1)THEN DD=P(IND) IND=IND+1 ELSE DD=DDM(J) ENDIF IF(B15(J).EQ.1)THEN CONC=P(IND) IND=IND+1 ELSE CONC=CONCM(J) ENDIF DO J=1,ACQ IF(B4(J).EQ.1)THEN TAUM=P(IND) IND=IND+1 ENDIF IF(B4(J).EQ.2)THEN TM1(1)=P(IND) TM2(1)=P(IND+1) IND=IND+2 ENDIF IF(B4(J).EQ.0) TAUM=TAUMM(J,1) IF(B5(J).EQ.1)THEN AMOLFRA=P(IND) IND=IND+1 ELSE AMOLFRA=AMOLFRAM(J) ENDIF IF(B6(J).EQ.1)THEN RK=P(IND) IND=IND+1 ELSE RK=RKM(J) ENDIF IF(B16(J).EQ.1)THEN ACONT=P(IND) IND=IND+1 ELSE ACONT=ACONTM(J) ENDIF IF(B4(J).EQ.2)TAUM=TM1(1)*EXP(TM2(1)/TEMP(K)) IF(B4(J).EQ.3)TAUM=TAUMM(J,K) C******************************************************************************* ADD=0 DO IJK=1,K-1 ADD=ADD+NPT(IJK) END DO IPLUS=I+ADD BZ=Z(IPLUS)*1000000. RK1=10.*2.46502D-52/(2.6752D8)**2*GAMMAI**2*(RK*1.E-10)**(-6) IF (AMOLFRA.EQ.0.)RK1=10./(SPIN*(SPIN+1.)*2./15.)/TAUS0*1.E-9 IF(RK.EQ.0) RK1=0. CONA=(ACONT*6.28*1.E6*1.0546E-34) C READ CALCULATED CONTRIBUTIONS TO TMUNO TMUNO=TMAT(IPLUS,J) TMUNOCONT=TMATCONT(IPLUS,J) TMUNOCROSS=TMATCROSS(IPLUS,J) C CALCULATION OF TMUNO TMUNO=TMUNO*RK1+CONA*CONA*TMUNOCONT+SQRT(RK1)*CONA*TMUNOCROSS TPUNO1=1./(1./TMUNO+TAUM)*AMOLFRA IF(AMOLFRA.EQ.0.)TPUNO1=1./(1./TMUNO+TAUM) TPUNO(IPLUS)=TPUNO1/CONC*0.001 TPUNOTOT=TPUNOTOT+TPUNO(IPLUS) END DO 56 CONTINUE C CALCULATION OF OUTER-SPHERE CONTRIBUTION TERM=0 IF(DD.NE.0)THEN V=BZ*2.*3.14 TAUD=D**2/DD CZ=SQRT(2*V*TAUD) ZZ=SQRT(2*V*TAUD*658.) GEI=(1.+5.*CZ/8.+CZ**2/8.)/(1.+CZ+CZ**2/2.+CZ**3/6.+4.* & CZ**4/81.+CZ**5/81.+CZ**6/648.) GES=(1.+5.*ZZ/8.+ZZ**2/8.)/(1.+ZZ+ZZ**2/2.+ZZ**3/6.+4.* & ZZ**4/81.+ZZ**5/81.+ZZ**6/648.) PRIMO=32.*3.14/405.*(2.6752E4*1.7608E7*1.0546E-27)**2*6.022E20 SECONDO=SPIN*(SPIN+1)*CONC/(D*DD) TERZO=(3.*GEI+7.*GES) TERM=(PRIMO*SECONDO*TERZO) ENDIF TPUNOTOT=TPUNOTOT+TERM C DIFFERENCES BETWEEN EXPERIMENTAL AND FITTED VALUES TPUNO(IPLUS)=TPUNOTOT FB=((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FB FBW=SQRT((PP(IPLUS)-TPUNO(IPLUS))**2)/PP(IPLUS)+FBW !IF (STEPGAMMA.NE.1) WRITE(6,'(2X,2(E10.4))')Z(IPLUS),TPUNO(IPLUS) IF(INDEXSTAMPA.EQ.1) WRITE(44,'(2X,2(E10.4))') Z(IPLUS), & TPUNO(IPLUS) END DO 223 CONTINUE FUNC=FBW RETURN END FUNCTION E(BZ,BETA,THETA,TAUC,NMX,PHI,GAMMA) IMPLICIT REAL*8(A-H,O-Z) COMMON /A3/ T11,T12,T13 COMMON /T1T2/ IREL COMMON /ECOM/ ECONT,ECROSS COMMON /STEPGAMMA/ STEPGAMMA OMI=BZ*6.2831853 CALL DIAG(BETA,GAMMA,BZ,NMX) IF(IREL.EQ.1.AND.STEPGAMMA.GT.1)CALL TUNO(BETA,OMI,THETA, & TAUC,NMX,PHI,GAMMA) IF(IREL.EQ.1.AND.STEPGAMMA.EQ.1)CALL TUNOISO(BETA,OMI,THETA, & TAUC,NMX) IF(IREL.EQ.2)CALL TDUE(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA) E=T11*SIN(BETA) ECONT=T12*SIN(BETA) ECROSS=T13*SIN(BETA) RETURN END SUBROUTINE DIAG(BETA,GAMMA,BZ,NMX) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(PI2 = 6.2831853, VL = 2.9979D+10) COMMON /RK10/ SPIN, SI COMMON /T1T2/ IREL COMMON /INDEX/ INDEX COMMON /STEPGAMMA/ STEPGAMMA COMMON /TAUDELTA/ TAUDELTA COMMON /TAUE/ TAUE COMMON /CONTAT/ ACONT COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4 COMMON /GTENSOR/ GX,GY,GZ COMMON /TAU/ TAUR,TAUV,TAUM COMMON /TAU1/ TAUS0 COMMON /AOLD/ OMOLD(10000),COLD(10000,4) COMPLEX*16 C(100,100,19) COMMON /A1/ OM(1000,1000),C PARAMETER(MBRANC=90) COMPLEX*16 SZ(MBRANC,MBRANC),SP(MBRANC,MBRANC),SM(MBRANC,MBRANC) DIMENSION CO1(MBRANC),CO2(MBRANC), CO3(MBRANC) COMPLEX*16 TZ,TP,TM DIMENSION WR( MBRANC ) COMPLEX*16 AR(MBRANC,MBRANC),ZR(MBRANC,MBRANC ) DIMENSION ARR(MBRANC,MBRANC),ARI(MBRANC,MBRANC) DIMENSION ZRR(MBRANC,MBRANC),ZRI(MBRANC,MBRANC) DIMENSION WK1(MBRANC),WK2(MBRANC),WK3(MBRANC) CT=COS(BETA) ST=SIN(BETA) C CALCULATION OF CORRELATION TIME WI=2*3.1416*BZ WS=658.2*WI IF(TAUV.EQ.0.AND.TAUR.EQ.0)THEN TAUC=TAUS0 IF(TAUM.EQ.0)THEN TAUE=TAUS0 ELSE TAUE=1./(1./TAUS0+1./TAUM) ENDIF ELSE STI=WS**2*TAUV**2 IF (TAUDELTA.EQ.2)THEN RTAUS=2.*(TAUS0*VL*PI2)**2*(4.*SPIN*(SPIN+1)-3)/50.* & (TAUV/(1+STI)+4*TAUV/(1+4*STI)) ELSE RTAUS=(0.2/TAUS0)*(1./(1.+STI)+4./(1.+4.*STI)) ENDIF IF(TAUR.EQ.0)THEN RTAUC=RTAUS ELSE RTAUC=RTAUS+1./TAUR ENDIF IF(TAUM.NE.0)THEN RTAUC=RTAUC+1./TAUM ENDIF TAUC=1./RTAUC IF (ACONT.NE.0)THEN IF (TAUM.EQ.0)THEN RTAUE=RTAUS ELSE RTAUE=RTAUS+1./TAUM ENDIF TAUE=1./RTAUE ENDIF ENDIF IF (STEPGAMMA.GT.1)THEN COEFFH=-1. ELSE COEFFH=1. ENDIF IF(ACONIND.EQ.1.)GO TO 456 IF (DPARATOT.EQ.0..AND.EPARATOT.EQ.0..AND.SPIN.EQ.0.5.AND. & GX.EQ.GY.AND.GX.EQ.GZ.AND.IREL.EQ.1)THEN C MATRIX OF ENERGY FOR HYPERFINE COUPLING X=BZ*3.1415927*658.2 ZC=X*CT ZS=X*ST DO 200 I=1,(2.*SI+1.)*2. DO 200 J=1,(2.*SI+1.)*2. 200 AR(I,J)=0. SSI = SI*(SI + 1.) DO I = 1, (2.*SI+1.)*2., 2 COEF = SI - (I - 1)/2 AR(I,I) = ZC*GZ + (SI-I/2)*AZ/2. AR(I+1,I+1) = -(ZC*GZ + (SI-I/2)*AZ/2.) AR(I,I+1) = COEFFH*ZS*GY AR(I+1,I) = COEFFH*ZS*GY AR(I+1,I+2) = 0.5*(AX+AY)/2.*SQRT(SSI-(COEF-1.)*COEF) AR(I+2,I+1) = 0.5*(AX+AY)/2.*SQRT(SSI-(COEF-1.)*COEF) END DO IF (INDEX.EQ.1)THEN WRITE(6,*)'DIM. MATRIX', NMX OPEN(UNIT=17,FILE='MAT') DO I=1,(2.*SI+1)*(2.*SPIN+1) DO J=1,(2.*SI+1)*(2.*SPIN+1) WRITE(17,*)AR(I,J) END DO WRITE(17,*)' ' END DO CLOSE(17) ENDIF INDEX=INDEX+1 C DIAGONALISATION OF THE MATRIX OF ENERGY DO 45 I=1,NMX DO 45 J=1,NMX ARR(I,J)=REAL(AR(I,J)) ARI(I,J)=IMAG(AR(I,J)) 45 CONTINUE CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC $ ,WK1,WK2,WK3,0) DO 46 I=1,NMX DO 46 J=1,NMX ZR(I,J)=CMPLX(ZRR(I,J),ZRI(I,J)) 46 CONTINUE I=1 OM(1,1)=0. OMOLD(1)=0. DO 700 K=1,NMX DO 700 L=1,NMX IF (K.EQ.L)GO TO 700 I=I+1 OM(K,L)=WR(K)-WR(L) C DIFFERENCES IN ENERGY LEVELS OMOLD(I)=WR(K)-WR(L) 700 CONTINUE C CALCULATION OF CORRELATION FUNCTIONS DO 400 K=1,NMX DO 400 L=1,NMX TZ=0 DO 1500 J=1,NMX TZ=-((-1.)**J)*ZR(J,K)*CONJG(ZR(J,L))+TZ 1500 CONTINUE SZ(K,L)=TZ/2. TP=0 DO J=1,NMX,2 TP=ZR(J,K)*CONJG(ZR(J+1,L))+TP END DO SP(K,L)=TP TM=0 DO J=2,NMX,2 TM=ZR(J,K)*CONJG(ZR(J-1,L))+TM END DO SM(K,L)=TM 400 CONTINUE GO TO 567 ENDIF IF (APARTOT.EQ.0.AND.APERTOT.EQ.0.AND.APERTOT2.EQ.0. & AND.EPARATOT.EQ.0.AND.GX.EQ.GY.AND.GX.EQ.GZ.AND.IREL.EQ.1)THEN IF (INDEX.EQ.1)THEN WRITE(6,*)'DIM. MATRIX', NMX ENDIF INDEX=INDEX+1 C MATRIX OF ENERGY IN ZERO FIELD SPLITTING X=BZ*2*3.1415927*658.2 ZC=X*CT ZS=X*ST DO 5200 I=1,NMX DO 5200 J=1,NMX 5200 AR(I,J)=0. S = FLOAT(NMX - 1)/2. SS = S*(S + 1.) DO I = 1, NMX COEF = S - DFLOAT(I - 1) AR(I,I) = COEF*ZC*GZ + DPARA*(COEF**2 - SS/3.) IF(I.LT.NMX) THEN AR(I,I+1) = COEFFH*0.5*ZS*GY*SQRT(SS-(COEF-1.)*COEF) AR(I+1,I) = AR(I,I+1) END IF END DO DO 145 I=1,NMX DO 145 J=1,NMX ARR(I,J)=REAL(AR(I,J)) ARI(I,J)=IMAG(AR(I,J)) 145 CONTINUE CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC $ ,WK1,WK2,WK3,0) DO 146 I=1,NMX DO 146 J=1,NMX ZR(I,J)=CMPLX(ZRR(I,J),ZRI(I,J)) 146 CONTINUE I=1 OM(1,1)=0. OMOLD(1)=0. DO 570 K=1,NMX DO 570 L=1,NMX IF (K.EQ.L)GO TO 570 I=I+1 OM(K,L)=WR(K)-WR(L) OMOLD(I)=WR(K)-WR(L) 570 CONTINUE C PER SPIN (SZ) DIVERSI DA 1/2 DO J=1,NMX CO1(J)=(NMX-(2*J-1))/2. END DO DO J=1,NMX-1 CO2(J)=SQRT(SS-CO1(J+1)*(CO1(J+1)+1.)) END DO DO J=2,NMX CO3(J)=SQRT(SS-CO1(J-1)*(CO1(J-1)-1.)) END DO DO 540 K=1, NMX DO 540 L=1, NMX TZ=0 DO 510 J=1,NMX TZ=CO1(J)*CONJG(ZR(J,K))*(ZR(J,L))+TZ 510 CONTINUE SZ(K,L)=TZ TP=0 DO 520 J=1,NMX-1 TP=CO2(J)*CONJG(ZR(J,K))*(ZR(J+1,L))+TP 520 CONTINUE SP(K,L)=TP TM=0 DO 530 J=2,NMX TM=CO3(J)*CONJG(ZR(J,K))*(ZR(J-1,L))+TM 530 CONTINUE SM(K,L)=TM 540 CONTINUE GO TO 567 ENDIF 456 CONTINUE C GENERAL MATRIX OF ENERGY X=BZ*3.1415927*658.2 ZC=X*CT ZS=X*ST DO I=1,(2.*SI+1)*(2.*SPIN+1)*2 DO J=1,(2.*SI+1)*(2.*SPIN+1)*2 AR(I,J)=0. END DO END DO SSI = SI*(SI + 1.) SS=SPIN*(SPIN+1.) ISMS=2.*SPIN+1 K=0 DO I=1,(2.*SI+1)*(2.*SPIN+1),2.*SPIN+1 DO J=0,2.*SPIN COEF = SI - K COEF2 = SPIN-J AR(I+J,I+J)=2.*(SPIN-J)*ZC*GZ + (SPIN-J)*AZ*(SI-K)+DPARA* $ (COEF2**2-SS/3.)+S4/120.*(35.*COEF2**4-30.*SS*COEF2**2+ $ 25.*COEF2**2-6.*SS+3.*SS**2) IF (J+1.LT.2.*SPIN+1) THEN AR(I+J,I+J+1)=CMPLX(COEFFH*ZS*GX*COS(GAMMA)*SQRT(SS- $ (COEF2-1.)*COEF2),-ZS*GY*SIN(GAMMA)*SQRT(SS-(COEF2-1.)*COEF2)) AR(I+J+1,I+J)=CMPLX(COEFFH*ZS*GX*COS(GAMMA)*SQRT(SS-(COEF2-1.)* $ COEF2),ZS*GY*SIN(GAMMA)*SQRT(SS-(COEF2-1.)*COEF2)) IF (J+2.LT.2.*SPIN+1) THEN AR(I+J,I+J+2)=EPARA* & SQRT(SS-COEF2*(COEF2-1.))* & SQRT(SS-(COEF2-1.)*(COEF2-2.))/2. AR(I+J+2,I+J)=AR(I+J,I+J+2) ENDIF IF (J+4.LT.2.*SPIN+1) THEN AR(I+J,I+J+4)=S4/48.* & SQRT(SS-COEF2*(COEF2-1.))* & SQRT(SS-(COEF2-1.)*(COEF2-2.))* & SQRT(SS-(COEF2-2.)*(COEF2-3.))* & SQRT(SS-(COEF2-3.)*(COEF2-4.)) AR(I+J+4,I+J)=AR(I+J,I+J+4) ENDIF IF (I+(2.*SPIN+1)+J.LE.(2.*SPIN+1)*(2.*SI+1))THEN AR(I+(ISMS)+J,(I+1+J))=.5*(AX+AY)/2.*SQRT(SSI-(COEF-1.)*COEF)* $ SQRT(SS-(COEF2-1.)*COEF2) AR((I+1+J),I+(ISMS)+J)= AR(I+(ISMS)+J,(I+1+J)) ENDIF IF (I+(2.*SPIN+1)+J+1.LE.(2.*SPIN+1)*(2.*SI+1))THEN AR(I+(ISMS)+J+1,(I+J))=.5*(AX-AY)/2.*SQRT(SSI-(COEF-1.)*COEF)* $ SQRT(SS-(COEF2-1.)*COEF2) AR((I+J),I+(ISMS)+J+1)= AR(I+(ISMS)+J+1,(I+J)) ENDIF ENDIF END DO K=K+1 END DO IF (INDEX.EQ.1)THEN WRITE(6,*)'DIM. MATRIX', NMX OPEN(UNIT=17,FILE='MAT') DO I=1,(2.*SI+1)*(2.*SPIN+1) DO J=1,(2.*SI+1)*(2.*SPIN+1) WRITE(17,*)AR(I,J) END DO WRITE(17,*)' ' END DO CLOSE(17) ENDIF INDEX=INDEX+1 C WR EIGENVALUES ZR EIGENVECTORS DO 245 I=1,NMX DO 245 J=1,NMX ARR(I,J)=REAL(AR(I,J)) ARI(I,J)=IMAG(AR(I,J)) 245 CONTINUE CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC $ ,WK1,WK2,WK3,0) DO 246 I=1,NMX DO 246 J=1,NMX ZR(I,J)=CMPLX(ZRR(I,J),ZRI(I,J)) 246 CONTINUE I=1 OM(1,1)=0. OMOLD(1)=0. DO 70 K=1,NMX DO 70 L=1,NMX IF (K.EQ.L)GO TO 70 I=I+1 OM(K,L)=WR(K)-WR(L) OMOLD(I)=WR(K)-WR(L) 70 CONTINUE J=0 DO I=1,(2.*SI+1)*(2.*SPIN+1) J=J+1 IF (J.GT.(2*SPIN+1)) J=1 CO1(I)=(2*SPIN+1-(2*J-1))/2. C CO1(I)=-((-1.)**J) END DO DO I=1,(2.*SI+1)*(2.*SPIN+1)-1 CO2(I)=SQRT(SPIN*(SPIN+1)-CO1(I+1)*(CO1(I+1)+1.)) END DO DO I=2,(2.*SI+1)*(2.*SPIN+1) CO3(I)=SQRT(SPIN*(SPIN+1)-CO1(I-1)*(CO1(I-1)-1.)) END DO IF (INDEX.EQ.2)THEN OPEN(UNIT=18,FILE='COE') DO I=1,NMX WRITE(18,*)CO1(I),CO2(I), CO3(I) END DO CLOSE(18) ENDIF DO 40 K=1,NMX DO 40 L=1, NMX TZ=CMPLX(0,0) DO 10 J=1,NMX TZ=CO1(J)*CONJG(ZR(J,K))*(ZR(J,L))+TZ 10 CONTINUE SZ(K,L)=TZ TP=CMPLX(0,0) DO I=1,(2.*SI+1)*(2.*SPIN+1),2.*SPIN+1 DO J=0,2.*SPIN-1 TP=CO2(I+J)*CONJG(ZR(I+J,K))*(ZR(I+J+1,L))+TP END DO END DO SP(K,L)=TP TM=CMPLX(0,0) DO I=1,(2.*SI+1)*(2.*SPIN+1),2.*SPIN+1 DO J=1,2.*SPIN TM=CO3(I+J)*CONJG(ZR(I+J,K))*(ZR(I+J-1,L))+TM END DO END DO SM(K,L)=TM 40 CONTINUE 567 CONTINUE IF (STEPGAMMA.GT.1)THEN I=1 C(1,1,11)=0. C(1,1,12)=0. C(1,1,13)=0. C(1,1,14)=0. C(1,1,15)=0. C(1,1,16)=0. C(1,1,17)=0. C(1,1,18)=0. C(1,1,19)=0. DO 150 K=1,NMX C AUTOCORRELATION FUNCTIONS C(1,1,11)=(GX**2/4.)*(SP(K,K)*SP(K,K) + SP(K,K)*SM(K,K) + & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K)) + & (GX*GY/2.)*(SP(K,K)*SP(K,K) - SM(K,K)*SM(K,K)) + & (GY**2/4.)*(SP(K,K)*SP(K,K) - SP(K,K)*SM(K,K) - & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K))+C(1,1,11) C(1,1,12)=(GX**2/4.)*(SP(K,K)*SP(K,K) + SP(K,K)*SM(K,K) + & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K)) & - (GX*GY/2.)*(SP(K,K)*SP(K,K) - SM(K,K)*SM(K,K)) + & (GY**2/4.)*(SP(K,K)*SP(K,K) - SP(K,K)*SM(K,K) - & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K))+C(1,1,12) C(1,1,13)=(GX**2/4.)*(SP(K,K)*SP(K,K) + SP(K,K)*SM(K,K) + & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K)) & - (GX*GY/2.)*(SP(K,K)*SM(K,K) - SM(K,K)*SP(K,K)) & - (GY**2/4.)*(SP(K,K)*SP(K,K) - SP(K,K)*SM(K,K) - & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K))+C(1,1,13) C(1,1,14)=(GX**2/4.)*(SP(K,K)*SP(K,K) + SP(K,K)*SM(K,K) + & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K)) & - (GX*GY/2.)*( - SP(K,K)*SM(K,K) + SM(K,K)*SP(K,K)) & - (GY**2/4.)*(SP(K,K)*SP(K,K) - SP(K,K)*SM(K,K) - & SM(K,K)*SP(K,K) + SM(K,K)*SM(K,K))+C(1,1,14) C(1,1,15)=(GX*GZ/2.)*(SP(K,K)*SZ(K,K) + SM(K,K)*SZ(K,K)) + & (GY*GZ/2.)*(SP(K,K)*SZ(K,K) - SM(K,K)*SZ(K,K))+C(1,1,15) C(1,1,16)=(GX*GZ/2.)*(SZ(K,K)*SP(K,K) + SZ(K,K)*SM(K,K)) + & (GY*GZ/2.)*(SZ(K,K)*SP(K,K) - SZ(K,K)*SM(K,K))+C(1,1,16) C(1,1,17)=(GX*GZ/2.)*(SP(K,K)*SZ(K,K) + SM(K,K)*SZ(K,K)) - & (GY*GZ/2.)*(SP(K,K)*SZ(K,K) - SM(K,K)*SZ(K,K))+C(1,1,17) C(1,1,18)=(GX*GZ/2.)*(SZ(K,K)*SP(K,K) + SZ(K,K)*SM(K,K)) - & (GY*GZ/2.)*(SZ(K,K)*SP(K,K) - SZ(K,K)*SM(K,K))+C(1,1,18) C(1,1,19)=GZ**2*SZ(K,K)*SZ(K,K)+C(1,1,19) 150 CONTINUE DO 100 K=1,NMX DO 100 L=1,NMX IF(K.EQ.L) GO TO 100 I=I+1 C(K,L,1)=(GX**2/4.)*(SP(K,L)*SP(L,K) + SP(K,L)*SM(L,K) + & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) + & (GX*GY/2.)*(SP(K,L)*SP(L,K) - SM(K,L)*SM(L,K)) + & (GY**2/4.)*(SP(K,L)*SP(L,K) - SP(K,L)*SM(L,K) - & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) C(K,L,2)=(GX**2/4.)*(SP(K,L)*SP(L,K) + SP(K,L)*SM(L,K) + & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) & - (GX*GY/2.)*(SP(K,L)*SP(L,K) - SM(K,L)*SM(L,K)) + & (GY**2/4.)*(SP(K,L)*SP(L,K) - SP(K,L)*SM(L,K) - & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) C(K,L,3)=(GX**2/4.)*(SP(K,L)*SP(L,K) + SP(K,L)*SM(L,K) + & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) & - (GX*GY/2.)*(SP(K,L)*SM(L,K) - SM(K,L)*SP(L,K)) & - (GY**2/4.)*(SP(K,L)*SP(L,K) - SP(K,L)*SM(L,K) - & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) C(K,L,4)=(GX**2/4.)*(SP(K,L)*SP(L,K) + SP(K,L)*SM(L,K) + & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) & - (GX*GY/2.)*( - SP(K,L)*SM(L,K) + SM(K,L)*SP(L,K)) & - (GY**2/4.)*(SP(K,L)*SP(L,K) - SP(K,L)*SM(L,K) - & SM(K,L)*SP(L,K) + SM(K,L)*SM(L,K)) C(K,L,5)=(GX*GZ/2.)*(SP(K,L)*SZ(L,K) + SM(K,L)*SZ(L,K)) + & (GY*GZ/2.)*(SP(K,L)*SZ(L,K) - SM(K,L)*SZ(L,K)) C(K,L,6)=(GX*GZ/2.)*(SZ(K,L)*SP(L,K) + SZ(K,L)*SM(L,K)) + & (GY*GZ/2.)*(SZ(K,L)*SP(L,K) - SZ(K,L)*SM(L,K)) C(K,L,7)=(GX*GZ/2.)*(SP(K,L)*SZ(L,K) + SM(K,L)*SZ(L,K)) - & (GY*GZ/2.)*(SP(K,L)*SZ(L,K) - SM(K,L)*SZ(L,K)) C(K,L,8)=(GX*GZ/2.)*(SZ(K,L)*SP(L,K) + SZ(K,L)*SM(L,K)) - & (GY*GZ/2.)*(SZ(K,L)*SP(L,K) - SZ(K,L)*SM(L,K)) C(K,L,9)=GZ**2*SZ(K,L)*SZ(L,K) 100 CONTINUE ELSE I=1 COLD(1,1)=0 COLD(1,2)=0 COLD(1,3)=0 COLD(1,4)=0 c gx=gy=gz=g/2.003 DO 1501 K=1,NMX C COEFFICIENTS FOR OMEGA_I COLD(1,1)=gz*gz*SP(K,K)**2+COLD(1,1) COLD(1,2)=gz*gz*SP(K,K)*SZ(K,K)+COLD(1,2) COLD(1,3)=gz*gz*SP(K,K)*SM(K,K)+COLD(1,3) COLD(1,4)=gz*gz*SZ(K,K)**2+COLD(1,4) 1501 CONTINUE DO 1001 K=1,NMX DO 1001 L=1,NMX IF(K.EQ.L) GO TO 1001 I=I+1 COLD(I,1)=gz*gz*SP(K,L)*SP(L,K) COLD(I,2)=gz*gz*SP(K,L)*SZ(L,K) COLD(I,3)=gz*gz*SP(K,L)*SM(L,K) COLD(I,4)=gz*gz*SZ(K,L)*SZ(L,K) 1001 CONTINUE ENDIF RETURN END SUBROUTINE POWELL(P,XI,N,NP,FTOL,ITER,FRET,NMX) C PERFORMES THE FITTING PROCEDURE IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=20,ITMAX=2000) DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX) COMMON /NMXCOM/ NMX1 COMMON /NVNP/ NV NMX1=NMX NV=NP CALL XISTEP(NV,XI,P) CALL FUNCZFS(P,FUNC,NMX,NP) FRET=FUNC DO 11 J=1,N PT(J)=P(J) 11 CONTINUE ITER=0 1 ITER=ITER+1 FP=FRET IBIG=0 DEL=0. DO 13 I=1,N DO 12 J=1,N XIT(J)=XI(J,I) 12 CONTINUE FPTT=FRET CALL LINMIN(P,XIT,N,FRET) IF(ABS(FPTT-FRET).GT.DEL)THEN DEL=ABS(FPTT-FRET) IBIG=I ENDIF 13 CONTINUE IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN IF(ITER.EQ.ITMAX) PAUSE 'POWELL EXCEEDING MAXIMUM ITERATIONS.' DO 14 J=1,N PTT(J)=2.*P(J)-PT(J) XIT(J)=P(J)-PT(J) PT(J)=P(J) 14 CONTINUE CALL FUNCZFS(PTT,FUNC,NMX,NP) FPTT=FUNC IF(FPTT.GE.FP)GO TO 1 T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2 IF(T.GE.0.)GO TO 1 CALL LINMIN(P,XIT,N,FRET) DO 15 J=1,N XI(J,IBIG)=XIT(J) 15 CONTINUE GO TO 1 END SUBROUTINE POWELLINT(P,XI,N,NP,FTOL,ITER,FRET,NMX) C PERFORMES THE INTERNAL FITTING PROCEDURE IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=20,ITMAX=2000) DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX) COMMON /NMXCOM/ NMX1 COMMON /NVNP2/ NV NMX1=NMX NV=NP CALL XISTEP2(NV,XI,P) CALL FUNCINT(P,FUNC,NMX,NP) FRET=FUNC DO 11 J=1,N PT(J)=P(J) 11 CONTINUE ITER=0 1 ITER=ITER+1 FP=FRET IBIG=0 DEL=0. DO 13 I=1,N DO 12 J=1,N XIT(J)=XI(J,I) 12 CONTINUE ! FPTT=FRET CALL LINMIN2(P,XIT,N,FRET) ! IF(ABS(FPTT-FRET).GT.DEL)THEN IF(ABS(FP-FRET).GT.DEL)THEN ! DEL=ABS(FPTT-FRET) DEL=ABS(FP-FRET) IBIG=I ENDIF 13 CONTINUE IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN IF(ITER.EQ.ITMAX) RETURN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF(ITER.EQ.ITMAX) PAUSE 'POWELL EXCEEDING MAXIMUM ITERATIONS.' DO 14 J=1,N PTT(J)=2.*P(J)-PT(J) XIT(J)=P(J)-PT(J) PT(J)=P(J) 14 CONTINUE CALL FUNCINT(PTT,FUNC,NMX,NP) FPTT=FUNC IF(FPTT.GE.FP)GO TO 1 T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2 IF(T.GE.0.)GO TO 1 CALL LINMIN2(P,XIT,N,FRET) DO 15 J=1,N XI(J,IBIG)=XIT(J) 15 CONTINUE GO TO 1 END SUBROUTINE LINMIN(P,XI,N,FRET) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=50,TOL=1.E-4) EXTERNAL F1DIM DIMENSION P(N),XI(N) COMMON /F1COM/ PCOM(NMAX),XICOM(NMAX),NCOM NCOM=N DO 11 J=1,N PCOM(J)=P(J) XICOM(J)=XI(J) 11 CONTINUE AX=0. XX=1 ! BX=2 CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM) FRET=BRENT(AX,XX,BX,F1DIM,TOL,XMIN) DO 12 J=1,N XI(J)=XMIN*XI(J) P(J)=P(J)+XI(J) 12 CONTINUE RETURN END SUBROUTINE LINMIN2(P,XI,N,FRET) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=50,TOL=1.E-4) EXTERNAL F2DIM DIMENSION P(N),XI(N) COMMON /F2COM/ PCOM(NMAX),XICOM(NMAX),NCOM NCOM=N DO 11 J=1,N PCOM(J)=P(J) XICOM(J)=XI(J) 11 CONTINUE AX=0. XX=1 ! BX=2 CALL MNBRAK2(AX,XX,BX,FA,FX,FB,F2DIM) FRET=BRENT2(AX,XX,BX,F2DIM,TOL,XMIN) DO 12 J=1,N XI(J)=XMIN*XI(J) P(J)=P(J)+XI(J) 12 CONTINUE RETURN END FUNCTION F1DIM(X) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=50) COMMON /F1COM/ PCOM(NMAX),XICOM(NMAX),NCOM COMMON /NMXCOM/ NMX1 COMMON /NVNP/ NV DIMENSION XT(NMAX) DO 11 J=1,NCOM XT(J)=PCOM(J)+X*XICOM(J) 11 CONTINUE CALL FUNCZFS(XT,FUNC,NMX1,NV) F1DIM=FUNC RETURN END FUNCTION F2DIM(X) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=50) COMMON /F2COM/ PCOM(NMAX),XICOM(NMAX),NCOM COMMON /NMXCOM/ NMX1 COMMON /NVNP2/ NV DIMENSION XT(NMAX) DO 11 J=1,NCOM XT(J)=PCOM(J)+X*XICOM(J) 11 CONTINUE CALL FUNCINT(XT,FUNC,NMX1,NV) F2DIM=FUNC RETURN END SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,F1DIM) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20) FA=F1DIM(AX) FB=F1DIM(BX) IF(FB.GT.FA)THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM ENDIF CX=BX+GOLD*(BX-AX) FC=F1DIM(CX) 1 IF(FB.GE.FC)THEN R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R)) ULIM=BX+GLIMIT*(CX-BX) IF((BX-U)*(U-CX).GT.0.)THEN FU=F1DIM(U) IF(FU.LT.FC)THEN AX=BX FA=FB BX=U FB=FU GO TO 1 ELSE IF(FU.GT.FB)THEN CX=U FC=FU GO TO 1 ENDIF U=CX+GOLD*(CX-BX) FU=F1DIM(U) ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN FU=F1DIM(U) IF(FU.LT.FC)THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU FU=F1DIM(U) ENDIF ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN U=ULIM FU=F1DIM(U) ELSE U=CX+GOLD*(CX-BX) FU=F1DIM(U) ENDIF AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GO TO 1 ENDIF RETURN END SUBROUTINE MNBRAK2(AX,BX,CX,FA,FB,FC,F2DIM) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20) FA=F2DIM(AX) FB=F2DIM(BX) IF(FB.GT.FA)THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM ENDIF CX=BX+GOLD*(BX-AX) FC=F2DIM(CX) 1 IF(FB.GE.FC)THEN R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R)) ULIM=BX+GLIMIT*(CX-BX) IF((BX-U)*(U-CX).GT.0.)THEN FU=F2DIM(U) IF(FU.LT.FC)THEN AX=BX FA=FB BX=U FB=FU GO TO 1 ELSE IF(FU.GT.FB)THEN CX=U FC=FU GO TO 1 ENDIF U=CX+GOLD*(CX-BX) FU=F2DIM(U) ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN FU=F2DIM(U) IF(FU.LT.FC)THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU FU=F2DIM(U) ENDIF ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN U=ULIM FU=F2DIM(U) ELSE U=CX+GOLD*(CX-BX) FU=F2DIM(U) ENDIF AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GO TO 1 ENDIF RETURN END FUNCTION BRENT(AX,BX,CX,F1DIM,TOL,XMIN) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10) A=MIN(AX,CX) B=MAX(AX,CX) V=BX W=V X=V E=0. FX=F1DIM(X) FV=FX FW=FX DO 11 ITER=1,ITMAX XM=0.5*(A+B) TOL1=TOL*ABS(X)+ZEPS TOL2=2.*TOL1 IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3 IF(ABS(E).GT.TOL1) THEN R=(X-W)*(FX-FV) Q=(X-V)*(FX-FW) P=(X-V)*Q-(X-W)*R Q=2.*(Q-R) IF(Q.GT.0.) P=-P Q=ABS(Q) ETEMP=E E=D IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR. * P.GE.Q*(B-X)) GOTO 1 D=P/Q U=X+D IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X) GOTO 2 ENDIF 1 IF(X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF D=CGOLD*E 2 IF(ABS(D).GE.TOL1) THEN U=X+D ELSE U=X+SIGN(TOL1,D) ENDIF FU=F1DIM(U) IF(FU.LE.FX) THEN IF(U.GE.X) THEN A=X ELSE B=X ENDIF V=W FV=FW W=X FW=FX X=U FX=FU ELSE IF(U.LT.X) THEN A=U ELSE B=U ENDIF IF(FU.LE.FW .OR. W.EQ.X) THEN V=W FV=FW W=U FW=FU ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN V=U FV=FU ENDIF ENDIF 11 CONTINUE ! PAUSE 'BRENT EXCEED MAXIMUM ITERATIONS.' WRITE(6,*) 'BRENT EXCEED MAXIMUM ITERATIONS.' RETURN !!!!!!!!!!!!!!!!!!!!! 3 XMIN=X BRENT=FX RETURN END FUNCTION BRENT2(AX,BX,CX,F2DIM,TOL,XMIN) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10) A=MIN(AX,CX) B=MAX(AX,CX) V=BX W=V X=V E=0. FX=F2DIM(X) FV=FX FW=FX DO 11 ITER=1,ITMAX XM=0.5*(A+B) TOL1=TOL*ABS(X)+ZEPS TOL2=2.*TOL1 IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3 IF(ABS(E).GT.TOL1) THEN R=(X-W)*(FX-FV) Q=(X-V)*(FX-FW) P=(X-V)*Q-(X-W)*R Q=2.*(Q-R) IF(Q.GT.0.) P=-P Q=ABS(Q) ETEMP=E E=D IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR. * P.GE.Q*(B-X)) GOTO 1 D=P/Q U=X+D IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X) GOTO 2 ENDIF 1 IF(X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF D=CGOLD*E 2 IF(ABS(D).GE.TOL1) THEN U=X+D ELSE U=X+SIGN(TOL1,D) ENDIF FU=F2DIM(U) IF(FU.LE.FX) THEN IF(U.GE.X) THEN A=X ELSE B=X ENDIF V=W FV=FW W=X FW=FX X=U FX=FU ELSE IF(U.LT.X) THEN A=U ELSE B=U ENDIF IF(FU.LE.FW .OR. W.EQ.X) THEN V=W FV=FW W=U FW=FU ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN V=U FV=FU ENDIF ENDIF 11 CONTINUE ! PAUSE 'BRENT EXCEED MAXIMUM ITERATIONS.' WRITE(6,*) 'BRENT2 EXCEED MAXIMUM ITERATIONS.' RETURN !!!!!!!!!!!!!!!!!!!!! 3 XMIN=X BRENT2=FX RETURN END SUBROUTINE XISTEP(NV,XI,P) DIMENSION XI(NV,NV),P(NV) IMPLICIT REAL*8(A-H,O-Z) COMMON /ALFASTEP/ ALFASTEP DO 13 I=1,NV DO 13 J=1,NV 13 XI(I,J)=ALFASTEP*P(I) DO 14 I=1,NV 14 XI(I,I)=-ALFASTEP*P(I) RETURN END SUBROUTINE XISTEP2(NV,XI,P) DIMENSION XI(NV,NV),P(NV) IMPLICIT REAL*8(A-H,O-Z) COMMON /ALFASTEP/ ALFASTEP DO 13 I=1,NV DO 13 J=1,NV 13 XI(I,J)=ALFASTEP*P(I) DO 14 I=1,NV 14 XI(I,I)=-ALFASTEP*P(I) RETURN END SUBROUTINE F01BCF(N,TOL,Z,IZ,W,IW,D,E,C,S) C MARK 3 RELEASE. NAG COPYRIGHT 1972. C MARK 4 REVISED. C MARK 4.5 REVISED C MARK 5C REVISED C MARK 6B REVISED IER-125 IER-127 (JUL 1978) C MARK 11 REVISED. VECTORISATION (JAN 1984). C MARK 11.5(F77) REVISED. (SEPT 1985.) C C C TRECX2 C F01BCF REDUCES A COMPLEX HERMITIAN MATRIX TO REAL C TRIDIAGONAL FORM FROM WHICH THE EIGENVALUES AND EIGENVECTORS C CAN BE FOUND USING SUBROUTINE F02AYF,(CXTQL2). THE HERMITIAN C MATRIX A=A(1) IS REDUCED TO THE TRIDIAGONAL MATRIX A(N-1) BY C N-2 UNITARY TRANSFORMATIONS. THE HOUSEHOLDER REDUCTION ITSELF C DOES NOT GIVE A REAL TRIDIAGONAL MATRIX, THE OFF-DIAGONAL C ELEMENTS ARE COMPLEX. THEY ARE SUBSEQUENTLY MADE REAL BY A C DIAGONAL TRANSFORMATION. C APRIL 1ST. 1972 C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION TOL INTEGER IW, IZ, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION C(N), D(N), E(N), S(N), W(IW,N), Z(IZ,N) C .. LOCAL SCALARS .. DOUBLE PRECISION CO, F, FI, FR, G, GI, GR, H, HH, R, SI INTEGER I, II, J, K, L C .. EXTERNAL SUBROUTINES .. EXTERNAL F01BCY, F01BCZ C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, SQRT C .. EXECUTABLE STATEMENTS .. DO 20 I = 1, N D(I) = Z(N,I) E(I) = -W(N,I) 20 CONTINUE IF (N.EQ.1) GO TO 540 DO 360 II = 2, N I = N - II + 2 L = I - 2 G = 0.0D0 FR = D(I-1) FI = E(I-1) IF (L.EQ.0) GO TO 60 DO 40 K = 1, L G = G + D(K)*D(K) + E(K)*E(K) 40 CONTINUE 60 H = G + FR*FR + FI*FI C L IS NOW I-1 L = L + 1 IF (ABS(FR)+ABS(FI).NE.0.0D0) GO TO 80 R = 0.0D0 CO = 1.0D0 C(I) = 1.0D0 SI = 0.0D0 S(I) = 0.0D0 GO TO 140 80 IF (ABS(FR).LT.ABS(FI)) GO TO 100 R = ABS(FR)*SQRT(1.0D0+(FI/FR)**2) GO TO 120 100 R = ABS(FI)*SQRT(1.0D0+(FR/FI)**2) 120 SI = FI/R S(I) = -SI CO = FR/R C(I) = CO 140 IF (G.LE.TOL) GO TO 280 G = -SQRT(H) E(I) = G C E(I) HAS ITS FINAL REAL VALUE H = H - R*G C S*S + SR D(I-1) = (R-G)*CO E(I-1) = (R-G)*SI DO 160 J = 1, L Z(J,I) = D(J) W(J,I) = E(J) 160 CONTINUE CALL F01BCZ(Z,IZ,W,IW,L,D,E,C,S) C FORM P DO 180 J = 1, L C(J) = C(J)/H S(J) = S(J)/H 180 CONTINUE FR = 0.0D0 DO 200 J = 1, L FR = FR + C(J)*D(J) + S(J)*E(J) 200 CONTINUE C FORM K HH = FR/(H+H) C FORM Q DO 220 J = 1, L C(J) = C(J) - HH*D(J) S(J) = S(J) - HH*E(J) 220 CONTINUE C NOW FORM REDUCED A DO 260 J = 1, L FR = D(J) FI = E(J) GR = C(J) GI = S(J) DO 240 K = J, L Z(K,J) = (((Z(K,J)-GR*D(K))-GI*E(K))-FR*C(K)) - FI*S(K) W(K,J) = (((W(K,J)-GR*E(K))+GI*D(K))-FR*S(K)) + FI*C(K) 240 CONTINUE D(J) = Z(L,J) Z(I,J) = 0.0D0 E(J) = -W(L,J) W(I,J) = 0.0D0 W(J,J) = 0.0D0 260 CONTINUE GO TO 340 280 E(I) = R H = 0.0D0 DO 300 J = 1, L Z(J,I) = D(J) W(J,I) = E(J) 300 CONTINUE DO 320 J = 1, L Z(I,J) = 0.0D0 D(J) = Z(I-1,J) W(I,J) = 0.0D0 E(J) = -W(I-1,J) 320 CONTINUE 340 D(I) = H 360 CONTINUE C WE NOW FORM THE PRODUCT OF THE C HOUSEHOLDER MATRICES, OVERWRITING C ON Z AND W DO 500 I = 2, N L = I - 1 Z(N,L) = Z(L,L) Z(L,L) = 1.0D0 W(N,L) = E(L) W(L,L) = 0.0D0 H = D(I) IF (H.EQ.0.0D0) GO TO 460 DO 380 K = 1, L D(K) = 0.0D0 E(K) = 0.0D0 380 CONTINUE CALL F01BCY(Z,IZ,W,IW,L,L,Z(1,I),W(1,I),D,E) DO 400 K = 1, L D(K) = D(K)/H E(K) = -E(K)/H 400 CONTINUE DO 440 J = 1, L DO 420 K = 1, L Z(K,J) = Z(K,J) - Z(K,I)*D(J) + W(K,I)*E(J) W(K,J) = W(K,J) - Z(K,I)*E(J) - W(K,I)*D(J) 420 CONTINUE 440 CONTINUE 460 DO 480 J = 1, L Z(J,I) = 0.0D0 W(J,I) = 0.0D0 480 CONTINUE 500 CONTINUE W(N,N) = E(N) DO 520 I = 1, N D(I) = Z(N,I) Z(N,I) = 0.0D0 E(I) = W(N,I) W(N,I) = 0.0D0 520 CONTINUE 540 Z(N,N) = 1.0D0 W(N,N) = 0.0D0 E(1) = 0.0D0 C NOW WE MULTIPLY BY THE C COSTHETA + I SINTHETA COLUMN C FACTORS CO = 1.0D0 SI = 0.0D0 IF (N.EQ.1) RETURN DO 580 I = 2, N F = CO*C(I) - SI*S(I) SI = CO*S(I) + SI*C(I) CO = F DO 560 J = 1, N F = Z(J,I)*CO - W(J,I)*SI W(J,I) = Z(J,I)*SI + W(J,I)*CO Z(J,I) = F 560 CONTINUE 580 CONTINUE RETURN END SUBROUTINE F01BCY(AR,IAR,AI,IAI,M,N,BR,BI,CR,CI) C MARK 11 RELEASE. NAG COPYRIGHT 1983. C MARK 11.5(F77) REVISED. (SEPT 1985.) C C COMPUTES C = C + (A**H)*B (COMPLEX) WHERE C A IS RECTANGULAR M BY N. C C MUST BE DISTINCT FROM B. C C C .. SCALAR ARGUMENTS .. INTEGER IAI, IAR, M, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(M), BR(M), CI(N), CR(N) C .. LOCAL SCALARS .. DOUBLE PRECISION XI, XR INTEGER I, J C .. EXECUTABLE STATEMENTS .. DO 40 I = 1, N XR = CR(I) XI = CI(I) DO 20 J = 1, M XR = XR + AR(J,I)*BR(J) + AI(J,I)*BI(J) XI = XI + AR(J,I)*BI(J) - AI(J,I)*BR(J) 20 CONTINUE CR(I) = XR CI(I) = XI 40 CONTINUE RETURN END SUBROUTINE F01BCZ(AR,IAR,AI,IAI,N,BR,BI,CR,CI) C MARK 11 RELEASE. NAG COPYRIGHT 1983. C MARK 11.5(F77) REVISED. (SEPT 1985.) C C COMPUTES C = A*B (COMPLEX) WHERE C A IS A HERMITIAN N-BY-N MATRIX, C WHOSE LOWER TRIANGLE IS STORED IN A. C C MUST BE DISTINCT FROM B. C C C .. SCALAR ARGUMENTS .. INTEGER IAI, IAR, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(N), BR(N), CI(N), CR(N) C .. LOCAL SCALARS .. DOUBLE PRECISION YI, YR INTEGER I, IP1, J, NM1 C .. EXECUTABLE STATEMENTS .. DO 20 I = 1, N CR(I) = 0.0D0 CI(I) = 0.0D0 20 CONTINUE IF (N.EQ.1) GO TO 100 NM1 = N - 1 DO 80 I = 1, NM1 DO 40 J = I, N CR(J) = CR(J) + AR(J,I)*BR(I) - AI(J,I)*BI(I) CI(J) = CI(J) + AR(J,I)*BI(I) + AI(J,I)*BR(I) 40 CONTINUE YR = CR(I) YI = CI(I) IP1 = I + 1 DO 60 J = IP1, N YR = YR + AR(J,I)*BR(J) + AI(J,I)*BI(J) YI = YI + AR(J,I)*BI(J) - AI(J,I)*BR(J) 60 CONTINUE CR(I) = YR CI(I) = YI 80 CONTINUE 100 CR(N) = CR(N) + AR(N,N)*BR(N) - AI(N,N)*BI(N) CI(N) = CI(N) + AR(N,N)*BI(N) + AI(N,N)*BR(N) RETURN END SUBROUTINE F02AXF(AR,IAR,AI,IAI,N,WR,VR,IVR,VI,IVI,WK1,WK2,WK3, * IFAIL) C MARK 3 RELEASE. NAG COPYRIGHT 1972. C MARK 4.5 REVISED C MARK 9 REVISED. IER-327 (SEP 1981). C MARK 11.5(F77) REVISED. (SEPT 1985.) C MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988). C C EIGENVALUES AND EIGENVECTORS OF A COMPLEX HERMITIAN MATRIX C 1ST APRIL 1972 C C .. PARAMETERS .. CHARACTER*6 SRNAME PARAMETER (SRNAME='F02AXF') C .. SCALAR ARGUMENTS .. INTEGER IAI, IAR, IFAIL, IVI, IVR, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION AI(IAI,N), AR(IAR,N), VI(IVI,N), VR(IVR,N), * WK1(N), WK2(N), WK3(N), WR(N) C .. LOCAL SCALARS .. DOUBLE PRECISION A, B, MAX, SQ, SUM, XXXX INTEGER I, ISAVE, J C .. LOCAL ARRAYS .. CHARACTER*1 P01REC(1) C .. EXTERNAL FUNCTIONS .. DOUBLE PRECISION X02AJF, X02AKF INTEGER P01ABF EXTERNAL X02AJF, X02AKF, P01ABF C .. EXTERNAL SUBROUTINES .. EXTERNAL F01BCF, F02AYF C .. INTRINSIC FUNCTIONS .. INTRINSIC SQRT C .. EXECUTABLE STATEMENTS .. ISAVE = IFAIL DO 40 I = 1, N IF (AI(I,I).NE.0.0D0) GO TO 140 DO 20 J = 1, I VR(I,J) = AR(I,J) VI(I,J) = AI(I,J) 20 CONTINUE 40 CONTINUE CALL F01BCF(N,X02AKF()/X02AJF(),VR,IVR,VI,IVI,WR,WK1,WK2,WK3) IFAIL = 1 CALL F02AYF(N,X02AJF(),WR,WK1,VR,IVR,VI,IVI,IFAIL) IF (IFAIL.EQ.0) GO TO 60 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC) RETURN C NORMALISE 60 DO 120 I = 1, N SUM = 0.0D0 MAX = 0.0D0 DO 80 J = 1, N SQ = VR(J,I)*VR(J,I) + VI(J,I)*VI(J,I) SUM = SUM + SQ IF (SQ.LE.MAX) GO TO 80 MAX = SQ A = VR(J,I) B = VI(J,I) 80 CONTINUE IF (SUM.EQ.0.0D0) GO TO 120 SUM = 1.0D0/SQRT(SUM*MAX) DO 100 J = 1, N SQ = SUM*(VR(J,I)*A+VI(J,I)*B) VI(J,I) = SUM*(VI(J,I)*A-VR(J,I)*B) VR(J,I) = SQ 100 CONTINUE 120 CONTINUE RETURN 140 IFAIL = P01ABF(ISAVE,2,SRNAME,0,P01REC) RETURN END SUBROUTINE F02AYF(N,EPS,D,E,Z,IZ,W,IW,IFAIL) C MARK 3 RELEASE. NAG COPYRIGHT 1972. C MARK 4 REVISED. C MARK 4.5 REVISED C MARK 9 REVISED. IER-326 (SEP 1981). C MARK 11.5(F77) REVISED. (SEPT 1985.) C C CXTQL2 C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS OF A C HERMITIAN MATRIX, WHICH HAS BEEN REDUCED TO A REAL C TRIDIAGONAL MATRIX, T, GIVEN WITH ITS DIAGONAL ELEMENTS IN C THE ARRAY D(N) AND ITS SUB-DIAGONAL ELEMENTS IN THE LAST N C - 1 STORES OF THE ARRAY E(N), USING QL TRANSFORMATIONS. THE C EIGENVALUES ARE OVERWRITTEN ON THE DIAGONAL ELEMENTS IN THE C ARRAY D IN ASCENDING ORDER. THE REAL AND IMAGINARY PARTS OF C THE EIGENVECTORS ARE FORMED IN THE ARRAYS Z,W(N,N) C RESPECTIVELY, OVERWRITING THE ACCUMULATED TRANSFORMATIONS AS C SUPPLIED BY THE SUBROUTINE F01BCF. THE SUBROUTINE WILL FAIL C IF ALL EIGENVALUES TAKE MORE THAN 30*N ITERATIONS C 1ST APRIL 1972 C C .. PARAMETERS .. CHARACTER*6 SRNAME PARAMETER (SRNAME='F02AYF') C .. SCALAR ARGUMENTS .. DOUBLE PRECISION EPS INTEGER IFAIL, IW, IZ, N C .. ARRAY ARGUMENTS .. DOUBLE PRECISION D(N), E(N), W(IW,N), Z(IZ,N) C .. LOCAL SCALARS .. DOUBLE PRECISION B, C, F, G, H, P, R, S INTEGER I, I1, II, ISAVE, J, K, L, M, M1 C .. LOCAL ARRAYS .. CHARACTER*1 P01REC(1) C .. EXTERNAL FUNCTIONS .. INTEGER P01ABF EXTERNAL P01ABF C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, SQRT C .. EXECUTABLE STATEMENTS .. ISAVE = IFAIL IF (N.EQ.1) GO TO 40 DO 20 I = 2, N E(I-1) = E(I) 20 CONTINUE 40 E(N) = 0.0D0 B = 0.0D0 F = 0.0D0 J = 30*N DO 300 L = 1, N H = EPS*(ABS(D(L))+ABS(E(L))) IF (B.LT.H) B = H C LOOK FOR SMALL SUB-DIAG ELEMENT DO 60 M = L, N IF (ABS(E(M)).LE.B) GO TO 80 60 CONTINUE 80 IF (M.EQ.L) GO TO 280 100 IF (J.LE.0) GO TO 400 J = J - 1 C FORM SHIFT G = D(L) H = D(L+1) - G IF (ABS(H).GE.ABS(E(L))) GO TO 120 P = H*0.5D0/E(L) R = SQRT(P*P+1.0D0) H = P + R IF (P.LT.0.0D0) H = P - R D(L) = E(L)/H GO TO 140 120 P = 2.0D0*E(L)/H R = SQRT(P*P+1.0D0) D(L) = E(L)*P/(1.0D0+R) 140 H = G - D(L) I1 = L + 1 IF (I1.GT.N) GO TO 180 DO 160 I = I1, N D(I) = D(I) - H 160 CONTINUE 180 F = F + H C QL TRANSFORMATION P = D(M) C = 1.0D0 S = 0.0D0 M1 = M - 1 DO 260 II = L, M1 I = M1 - II + L G = C*E(I) H = C*P IF (ABS(P).LT.ABS(E(I))) GO TO 200 C = E(I)/P R = SQRT(C*C+1.0D0) E(I+1) = S*P*R S = C/R C = 1.0D0/R GO TO 220 200 C = P/E(I) R = SQRT(C*C+1.0D0) E(I+1) = S*E(I)*R S = 1.0D0/R C = C/R 220 P = C*D(I) - S*G D(I+1) = H + S*(C*G+S*D(I)) C FORM VECTOR DO 240 K = 1, N H = Z(K,I+1) Z(K,I+1) = S*Z(K,I) + C*H Z(K,I) = C*Z(K,I) - S*H H = W(K,I+1) W(K,I+1) = S*W(K,I) + C*H W(K,I) = C*W(K,I) - S*H 240 CONTINUE 260 CONTINUE E(L) = S*P D(L) = C*P IF (ABS(E(L)).GT.B) GO TO 100 280 D(L) = D(L) + F 300 CONTINUE C ORDER EIGEN VALUES ANS EIGENVECTORS DO 380 I = 1, N K = I P = D(I) I1 = I + 1 IF (I1.GT.N) GO TO 340 DO 320 J = I1, N IF (D(J).GE.P) GO TO 320 K = J P = D(J) 320 CONTINUE 340 IF (K.EQ.I) GO TO 380 D(K) = D(I) D(I) = P DO 360 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P P = W(J,I) W(J,I) = W(J,K) W(J,K) = P 360 CONTINUE 380 CONTINUE IFAIL = 0 RETURN 400 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC) RETURN END INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC) C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C MARK 13 REVISED. IER-621 (APR 1988). C MARK 13B REVISED. IER-668 (AUG 1988). C C P01ABF IS THE ERROR-HANDLING ROUTINE FOR THE NAG LIBRARY. C C P01ABF EITHER RETURNS THE VALUE OF IERROR THROUGH THE ROUTINE C NAME (SOFT FAILURE), OR TERMINATES EXECUTION OF THE PROGRAM C (HARD FAILURE). DIAGNOSTIC MESSAGES MAY BE OUTPUT. C C IF IERROR = 0 (SUCCESSFUL EXIT FROM THE CALLING ROUTINE), C THE VALUE 0 IS RETURNED THROUGH THE ROUTINE NAME, AND NO C MESSAGE IS OUTPUT C C IF IERROR IS NON-ZERO (ABNORMAL EXIT FROM THE CALLING ROUTINE), C THE ACTION TAKEN DEPENDS ON THE VALUE OF IFAIL. C C IFAIL = 1: SOFT FAILURE, SILENT EXIT (I.E. NO MESSAGES ARE C OUTPUT) C IFAIL = -1: SOFT FAILURE, NOISY EXIT (I.E. MESSAGES ARE OUTPUT) C IFAIL =-13: SOFT FAILURE, NOISY EXIT BUT STANDARD MESSAGES FROM C P01ABF ARE SUPPRESSED C IFAIL = 0: HARD FAILURE, NOISY EXIT C C FOR COMPATIBILITY WITH CERTAIN ROUTINES INCLUDED BEFORE MARK 12 C P01ABF ALSO ALLOWS AN ALTERNATIVE SPECIFICATION OF IFAIL IN WHICH C IT IS REGARDED AS A DECIMAL INTEGER WITH LEAST SIGNIFICANT DIGITS C CBA. THEN C C A = 0: HARD FAILURE A = 1: SOFT FAILURE C B = 0: SILENT EXIT B = 1: NOISY EXIT C C EXCEPT THAT HARD FAILURE NOW ALWAYS IMPLIES A NOISY EXIT. C C S.HAMMARLING, M.P.HOOPER AND J.J.DU CROZ, NAG CENTRAL OFFICE. C C .. SCALAR ARGUMENTS .. INTEGER IERROR, IFAIL, NREC CHARACTER*(*) SRNAME C .. ARRAY ARGUMENTS .. CHARACTER*(*) REC(*) C .. LOCAL SCALARS .. INTEGER I, NERR CHARACTER*72 MESS C .. EXTERNAL SUBROUTINES .. EXTERNAL P01ABZ, X04AAF, X04BAF C .. INTRINSIC FUNCTIONS .. INTRINSIC ABS, MOD C .. EXECUTABLE STATEMENTS .. IF (IERROR.NE.0) THEN C ABNORMAL EXIT FROM CALLING ROUTINE IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. IFAIL.EQ.-13 .OR. * (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN C NOISY EXIT CALL X04AAF(0,NERR) DO 20 I = 1, NREC CALL X04BAF(NERR,REC(I)) 20 CONTINUE IF (IFAIL.NE.-13) THEN WRITE (MESS,FMT=99999) SRNAME, IERROR CALL X04BAF(NERR,MESS) IF (ABS(MOD(IFAIL,10)).NE.1) THEN C HARD FAILURE CALL X04BAF(NERR, * ' ** NAG HARD FAILURE - EXECUTION TERMINATED' * ) CALL P01ABZ ELSE C SOFT FAILURE CALL X04BAF(NERR, * ' ** NAG SOFT FAILURE - CONTROL RETURNED') END IF END IF END IF END IF P01ABF = IERROR RETURN C 99999 FORMAT (' ** ABNORMAL EXIT FROM NAG LIBRARY ROUTINE ',A,': IFAIL', * ' =',I6) END SUBROUTINE P01ABZ C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C C TERMINATES EXECUTION WHEN A HARD FAILURE OCCURS. C C ******************** IMPLEMENTATION NOTE ******************** C THE FOLLOWING STOP STATEMENT MAY BE REPLACED BY A CALL TO AN C IMPLEMENTATION-DEPENDENT ROUTINE TO DISPLAY A MESSAGE AND/OR C TO ABORT THE PROGRAM. C ************************************************************* C .. EXECUTABLE STATEMENTS .. STOP END DOUBLE PRECISION FUNCTION X02AJF() C MARK 12 RELEASE. NAG COPYRIGHT 1986. C C RETURNS (1/2)*B**(1-P) IF ROUNDS IS .TRUE. C RETURNS B**(1-P) OTHERWISE C C .. EXECUTABLE STATEMENTS .. X02AJF = 0.11102230246251568000D-015 RETURN END DOUBLE PRECISION FUNCTION X02AKF() C MARK 12 RELEASE. NAG COPYRIGHT 1986. C C RETURNS B**(EMIN-1) (THE SMALLEST POSITIVE MODEL NUMBER) C C .. EXECUTABLE STATEMENTS .. X02AKF = 0.22250738585072014000D-307 RETURN END SUBROUTINE X04AAF(I,NERR) C MARK 7 RELEASE. NAG COPYRIGHT 1978 C MARK 7C REVISED IER-190 (MAY 1979) C MARK 11.5(F77) REVISED. (SEPT 1985.) C IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER C (STORED IN NERR1). C IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO C VALUE SPECIFIED BY NERR. C C .. SCALAR ARGUMENTS .. INTEGER I, NERR C .. LOCAL SCALARS .. INTEGER NERR1 C .. SAVE STATEMENT .. SAVE NERR1 C .. DATA STATEMENTS .. DATA NERR1/0/ C .. EXECUTABLE STATEMENTS .. IF (I.EQ.0) NERR = NERR1 IF (I.EQ.1) NERR1 = NERR RETURN END SUBROUTINE X04BAF(NOUT,REC) C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C C X04BAF WRITES THE CONTENTS OF REC TO THE UNIT DEFINED BY NOUT. C C TRAILING BLANKS ARE NOT OUTPUT, EXCEPT THAT IF REC IS ENTIRELY C BLANK, A SINGLE BLANK CHARACTER IS OUTPUT. C IF NOUT.LT.0, I.E. IF NOUT IS NOT A VALID FORTRAN UNIT IDENTIFIER, C THEN NO OUTPUT OCCURS. C C .. SCALAR ARGUMENTS .. INTEGER NOUT CHARACTER*(*) REC C .. LOCAL SCALARS .. INTEGER I C .. INTRINSIC FUNCTIONS .. INTRINSIC LEN C .. EXECUTABLE STATEMENTS .. IF (NOUT.GE.0) THEN C REMOVE TRAILING BLANKS DO 20 I = LEN(REC), 2, -1 IF (REC(I:I).NE.' ') GO TO 40 20 CONTINUE C WRITE RECORD TO EXTERNAL FILE 40 WRITE (NOUT,FMT=99999) REC(1:I) END IF RETURN C 99999 FORMAT (A) END SUBROUTINE GAUINT (BZ,TAUM,NMX) C PERFORMES THE INTEGRATION ON SPACE IMPLICIT REAL*8(A-H,O-Z) COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4 COMMON /TOT/ DPARATOT,EPARATOT,APARTOT,APERTOT,APERTOT2,ACONIND COMMON /GTENSOR/ GX,GY,GZ COMMON /STEPGAMMA/ STEPGAMMA COMMON /RK10/ SPIN, SI COMMON /TAU1/ TAUS0 COMMON /ECOM/ ECONT,ECROSS COMMON /T1T2/ IREL COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS DIMENSION ARRAY(20),T(20) DATA ARRAY/1.,1.,.555556,.888889,.555556,.347855,.652145,.652145,. 1347855,.236927,.478629,.568889,.478629,.236927,.171324,.360762,.46 27914,.467914,.360762,.171324/ DATAT/-.577350,.577350,-.774597,0.,.774597,-.861136,-.339981,.3399 181,.861136,-.906180,-.538469,0.,.538469,.906180,-.932470,-.661209, 2-.238619,.238619,.661209,.932470/ TMUNO=0. TMUNOCONT=0. TMUNOCROSS=0. TMUNOTOT=0. TMUNOTOTCONT=0. TMUNOTOTCROSS=0. GAMMA=0. IF(EPARA.NE.0.OR.GZ.NE.GX.OR.GZ.NE.GY.OR.AX.NE.AY)THEN STEPGAMMA=20. ELSE STEPGAMMA=1. ENDIF IF(IREL.NE.1)STEPGAMMA=20. IF(ACONIND.NE.0)STEPGAMMA=20. 55 A=0 H=1.5708/5. HP=H/2. DO 100 I=1,5 G=(2.*A+H)/2.0 A=A+H SP=0 SPCONT=0 SPCROSS=0 DO 50 J=10,14 BETA=HP*T(J)+G SP=SP+ARRAY(J)*E(BZ,BETA,THETA,TAUC,NMX,PHI,GAMMA) SPCONT=SPCONT+ARRAY(J)*ECONT SPCROSS=SPCROSS+ARRAY(J)*ECROSS 50 CONTINUE TMUNO=TMUNO+HP*SP TMUNOCONT=TMUNOCONT+HP*SPCONT TMUNOCROSS=TMUNOCROSS+HP*SPCROSS 100 CONTINUE TMUNOTOT=TMUNOTOT+TMUNO/STEPGAMMA TMUNOTOTCONT=TMUNOTOTCONT+TMUNOCONT/STEPGAMMA TMUNOTOTCROSS=TMUNOTOTCROSS+TMUNOCROSS/STEPGAMMA TMUNO=0. TMUNOCONT=0. TMUNOCROSS=0. GAMMA=GAMMA+6.28/STEPGAMMA IF(GAMMA.GE.6.28)THEN TMUNO=TMUNOTOT TMUNOCONT=TMUNOTOTCONT TMUNOCROSS=TMUNOTOTCROSS GOTO 56 ENDIF GOTO 55 56 CONTINUE RETURN END SUBROUTINE TUNO(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA) C FOUND CONTRIBUTIONS TO T1 TO BE INTEGRATED IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 C(100,100,19) COMMON /A1/ OM(1000,1000),C COMMON /RK10/ SPIN, SI COMMON /TAU1/ TAUS0 COMMON /TAUE/ TAUE COMMON /MOLFRAZ/ AMOLFRA COMMON /CONTAT/ ACONT COMPLEX*16 F0,F1,F2,FMU,FMD COMPLEX*16 A,B,CI,D,EI,F,G,H,AI COMPLEX*16 FF1,FF2,FF3,FF4,FF5,FF6,FF7,FF8,FF9 COMPLEX*16 GG1,GG2,GG3,GG4,GG5,GG6,GG7,GG8,GG9 COMPLEX*16 HH1,HH2,HH3,HH4,HH5,HH6,HH7,HH8,HH9 COMPLEX*16 CCF1,CCF2,CCF3,CCF4,CCF5,CCF6,CCF7,CCF8,CCF9 COMPLEX*16 CCG1,CCG2,CCG3,CCG4,CCG5,CCG6,CCG7,CCG8,CCG9 COMPLEX*16 CCH1,CCH2,CCH3,CCH4,CCH5,CCH6,CCH7,CCH8,CCH9 COMPLEX*16 AIA,AJ,ACF,ACG,ACH COMMON /A3/ T11,T12,T13 CT=COS(BETA) ST=SIN(BETA) C CONVERT DEG. ---> RAD (CA) CONVER = ATAN(1.0)/45.0 CA=COS(THETA* CONVER) SA=SIN(THETA* CONVER) CF=COS(PHI* CONVER) CF2=COS(2.*PHI* CONVER) SF=SIN(PHI* CONVER) SF2=SIN(2.*PHI* CONVER) C **************** RACAH'S NORMALIZED SPHERICAL HARMONICS F0=-(3.*CA**2-1.)/2. F1=SQRT(6.)/2.*SA*CA*CMPLX(CF,SF) FMU=-SQRT(6.)/2.*SA*CA*CMPLX(CF,-SF) F2=-SQRT(6.)/4.*SA**2*CMPLX(CF2,SF2) FMD=-SQRT(6.)/4.*SA**2*CMPLX(CF2,-SF2) C ELEMENTS OF THE ROTATION MATRIX A=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),SIN(2.*GAMMA)) B=(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA)) CI=(1.+CT)**2/4. D=-(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA)) EI=-ST**2/2. F=-(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA)) G=(1.-CT)**2/4. H=(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA)) AI=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),-SIN(2.*GAMMA)) C ************************ GA=1/TAUC AKL=0. DO 10 K=1,NMX DO 10 L=1,NMX IF(K.EQ.L) GO TO 10 c S+S+ FF1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.) & + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*(3./20.) & + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.)) & *C(K,L,1) c S-S- FF2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.) & + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*(3./20.) & + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.)) & *C(K,L,2) c S-S+ FF3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.) & /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.) & -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(K,L,3) c S+S- FF4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.) & /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.) & -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(K,L,4) c S+SZ FF5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.) & + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+ & G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(K,L,5) c SZS+ FF6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.) & + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+ & CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(K,L,6) c S-SZ FF7=-(AI*F0*FMU*(1./10.)*SQRT(1.5) + H*F0*F0* & SQRT(1./50.) + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.) & /SQRT(2.)+ EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+ & CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(K,L,7) c SZS- FF8=-(AI*F0*FMU*(1./10.)*SQRT(1.5) + F*F0*F0* & SQRT(1./50.)+ CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.) & /SQRT(2.)+ EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+ & G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(K,L,8) c SZSZ FF9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)* & CI*F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU* & SQRT(6./50.)+ (3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)* & AI*FMU*FMU)*C(K,L,9) c S+S+ GG1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.) & + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*(3./20.) & + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.)) & *C(L,K,1) c S-S- GG2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.) & + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*(3./20.) & + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.)) & *C(L,K,2) c S-S+ GG3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.) & /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.) & -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(L,K,3) c S+S- GG4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.) & /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.) & -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(L,K,4) c S+SZ GG5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.) & + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+ & G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(L,K,5) c SZS+ GG6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.) & + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+ & CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(L,K,6) c S-SZ GG7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.) & + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+ & CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(L,K,7) c SZS- GG8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.) & + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+ & G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(L,K,8) c SZSZ GG9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI*F1 & *FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT(6./50.)+ & (3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU) & *C(L,K,9) AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX) AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX) TMP=(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))* & 1.E+9/((OM(K,L)+OMI)**2+GA**2)+ & (REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))* & 1.E+9/((OM(L,K)-OMI)**2+GA**2) AKL=TMP+AKL 10 CONTINUE c S+S+ HH1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.) & + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*(3./20.) & + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.)) & *C(1,1,11) c S-S- HH2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.) & + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*(3./20.) & + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.)) & *C(1,1,12) c S-S+ HH3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.) & /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.) & -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(1,1,13) c S+S- HH4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.) & /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.) & -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(1,1,14) c S+SZ HH5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.) & + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+ & G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(1,1,15) c SZS+ HH6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.) & + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+ & CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(1,1,16) c S-SZ HH7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.) & + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+ & CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(1,1,17) c SZS- HH8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.) & + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+ & G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(1,1,18) c SZSZ HH9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI*F1 & *FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT(6./50.)+ & (3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU) & *C(1,1,19) ALL=2.*(REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX)) & *GA*1.E+9/(OMI**2+GA**2)+(IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7 & +HH8+HH9)/FLOAT(NMX)*OMI-IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7+ & HH8+HH9)/FLOAT(NMX)*OMI) & *1.E+9/(OMI**2+GA**2) T1SC=-ALL-AKL C DIPOLAR TERM T11=T1SC C CONTACT TERM IF (ACONT.NE.0)THEN IF(AMOLFRA.EQ.0)THEN RKCONT=1./(SPIN*(SPIN+1.)*2./3./ & (1.0546)**2*1.E34*1.E34* & (ACONT*6.28*1.0546E-34*1.E6)**2)/TAUS0 ELSE RKCONT=1. ENDIF C **************** F0=-(3.*CA**2-1.)/2. F1=SQRT(6.)/2.*SA*CA*CMPLX(CF,SF) FMU=-SQRT(6.)/2.*SA*CA*CMPLX(CF,-SF) F2=-SQRT(6.)/4.*SA**2*CMPLX(CF2,SF2) FMD=-SQRT(6.)/4.*SA**2*CMPLX(CF2,-SF2) C ************************ GA=1/TAUE GACR=1/TAUC T12FG=0. T13FG=0. DO 101 K=1,NMX DO 101 L=1,NMX IF(K.EQ.L) GO TO 101 c S+S+ FF1=(A/2.) & *C(K,L,1) c S-S- FF2=(AI/2) & *C(K,L,2) c S+S- FF3=-G/2. & *C(K,L,3) c S-S+ FF4=-CI/2 & *C(K,L,4) c SZS+ FF5=-B*SQRT(1./2.) & *C(K,L,5) c S+SZ FF6=-D*SQRT(1./2.) & *C(K,L,6) c SZS- FF7=H*SQRT(1./2.) & *C(K,L,7) c S-SZ FF8=F*SQRT(1./2.) & *C(K,L,8) c SZSZ FF9=EI & *C(K,L,9) c S+S+ GG1=(A/2.) & *C(L,K,1) c S-S- GG2=(AI/2) & *C(L,K,2) c S+S- GG3=-G/2. & *C(L,K,3) c S-S+ GG4=-CI/2 & *C(L,K,4) c SZS+ GG5=-B*SQRT(1./2.) & *C(L,K,5) c S+SZ GG6=-D*SQRT(1./2.) & *C(L,K,6) c SZS- GG7=H*SQRT(1./2.) & *C(L,K,7) c S-SZ GG8=F*SQRT(1./2.) & *C(L,K,8) c SZSZ GG9=EI & *C(L,K,9) C CROSSRELAXATION c S+S+ CCF1= (2*SQRT(1./40.)*F0*A + SQRT(3./40.)*FMU*(B+D) & + SQRT(3./20.)*FMD*(CI+G))*C(K,L,1) c S-S- CCF2= (SQRT(3./20.)*F2*(CI+G) + SQRT(3./40.)*F1*(H+F) & + 2*SQRT(1./40.)*F0*AI)*C(K,L,2) c S-S+ CCF3= (-2*SQRT(3./20.)*F2*A - SQRT(3./40.)*F1*(B+D) & -SQRT(1./40.)*F0*(CI+G))*C(K,L,3) c S+S- CCF4= (-SQRT(1./40.)*F0*(CI+G) - SQRT(3./40.)*FMU*(H+F) & -2*SQRT(3./20.)*FMD*AI)*C(K,L,4) c S+SZ CCF5= (-SQRT(1./20.)*F0*(B+D) - 2*SQRT(3./20.)*FMU*EI & -SQRT(3./10.)*FMD*(H+F))*C(K,L,5) c SZS+ CCF6= (2*SQRT(3./20.)*F1*A + SQRT(2./10.)*F0*(B+D) & +SQRT(3./20.)*FMU*(CI+G))*C(K,L,6) c S-SZ CCF7= (SQRT(3./10.)*F2*(B+D) + 2*SQRT(3./20.)*F1*EI & +SQRT(1./20.)*F0*(H+F))*C(K,L,7) c SZS- CCF8= (-SQRT(3./20.)*F1*(CI+G) - SQRT(2./10.)*F0*(H+F) & -2*SQRT(3./20.)*FMU*AI)*C(K,L,8) c SZSZ CCF9= (-SQRT(3./10.)*F1*(B+D) - 2*SQRT(2./5.)*F0*EI & -SQRT(3./10.)*FMU*(H+F))*C(K,L,9) c S+S+ CCG1= (2*SQRT(1./40.)*F0*A + SQRT(3./40.)*FMU*(B+D) & + SQRT(3./20.)*FMD*(CI+G))*C(L,K,1) c S-S- CCG2= (SQRT(3./20.)*F2*(CI+G) + SQRT(3./40.)*F1*(H+F) & + 2*SQRT(1./40.)*F0*AI)*C(L,K,2) c S-S+ CCG3= (-2*SQRT(3./20.)*F2*A - SQRT(3./40.)*F1*(B+D) & -SQRT(1./40.)*F0*(CI+G))*C(L,K,3) c S+S- CCG4= (-SQRT(1./40.)*F0*(CI+G) - SQRT(3./40.)*FMU*(H+F) & -2*SQRT(3./20.)*FMD*AI)*C(L,K,4) c S+SZ CCG5= (-SQRT(1./20.)*F0*(B+D) - 2*SQRT(3./20.)*FMU*EI & -SQRT(3./10.)*FMD*(H+F))*C(L,K,5) c SZS+ CCG6= (2*SQRT(3./20.)*F1*A + SQRT(2./10.)*F0*(B+D) & +SQRT(3./20.)*FMU*(CI+G))*C(L,K,6) c S-SZ CCG7= (SQRT(3./10.)*F2*(B+D) + 2*SQRT(3./20.)*F1*EI & +SQRT(1./20.)*F0*(H+F))*C(L,K,7) c SZS- CCG8= (-SQRT(3./20.)*F1*(CI+G) - SQRT(2./10.)*F0*(H+F) & -2*SQRT(3./20.)*FMU*AI)*C(L,K,8) c SZSZ CCG9= (-SQRT(3./10.)*F1*(B+D) - 2*SQRT(2./5.)*F0*EI & -SQRT(3./10.)*FMU*(H+F))*C(L,K,9) AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX) AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX) TMP=(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))* & 1.E+34/((OM(K,L)+OMI)**2+GA**2)+ & (REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))* & 1.E+34/((OM(L,K)-OMI)**2+GA**2) ACF=(CCF1+CCF2+CCF3+CCF4+ & CCF5+CCF6+CCF7+CCF8+CCF9)/FLOAT(NMX) ACG=(CCG1+CCG2+CCG3+CCG4+ & CCG5+CCG6+CCG7+CCG8+CCG9)/FLOAT(NMX) CRRFG=(REAL(ACF)*GACR+IMAG(ACF)*(OMI+OM(K,L)))* & 1.E+34*SQRT(1.E9)/((OM(K,L)+OMI)**2+GACR**2)+ & (REAL(ACG)*GACR+IMAG(ACG)*(-OMI+OM(L,K)))* & 1.E+34*SQRT(1.E9)/((OM(L,K)-OMI)**2+GACR**2) C CONTACT TERM T12FG=-TMP/(1.0546)**2*1.E34*RKCONT+T12FG T13FG=-CRRFG/1.0546+T13FG 101 CONTINUE c S+S+ CCH1= (2*SQRT(1./40.)*F0*A + SQRT(3./40.)*FMU*(B+D) & + SQRT(3./20.)*FMD*(CI+G))*C(1,1,11) c S-S- CCH2= (SQRT(3./20.)*F2*(CI+G) + SQRT(3./40.)*F1*(H+F) & + 2*SQRT(1./40.)*F0*AI)*C(1,1,12) c S-S+ CCH3= (-2*SQRT(3./20.)*F2*A - SQRT(3./40.)*F1*(B+D) & -SQRT(1./40.)*F0*(CI+G))*C(1,1,13) c S+S- CCH4= (-SQRT(1./40.)*F0*(CI+G) - SQRT(3./40.)*FMU*(H+F) & -2*SQRT(3./20.)*FMD*AI)*C(1,1,14) c S+SZ CCH5= (-SQRT(1./20.)*F0*(B+D) - 2*SQRT(3./20.)*FMU*EI & -SQRT(3./10.)*FMD*(H+F))*C(1,1,15) c SZS+ CCH6= (2*SQRT(3./20.)*F1*A + SQRT(2./10.)*F0*(B+D) & +SQRT(3./20.)*FMU*(CI+G))*C(1,1,16) c S-SZ CCH7= (SQRT(3./10.)*F2*(B+D) + 2*SQRT(3./20.)*F1*EI & +SQRT(1./20.)*F0*(H+F))*C(1,1,17) c SZS- CCH8= (-SQRT(3./20.)*F1*(CI+G) - SQRT(2./10.)*F0*(H+F) & -2*SQRT(3./20.)*FMU*AI)*C(1,1,18) c SZSZ CCH9= (-SQRT(3./10.)*F1*(B+D) - 2*SQRT(2./5.)*F0*EI & -SQRT(3./10.)*FMU*(H+F))*C(1,1,19) ACH=2*(CCH1+CCH2+CCH3+CCH4+ & CCH5+CCH6+CCH7+CCH8+CCH9)/FLOAT(NMX) CRRH=REAL(ACH)*GACR*1.E+34*SQRT(1.E9)/(OMI**2+GACR**2) T13FG=T13FG-CRRH/(1.0546) T12=T12FG T13=T13FG ENDIF RETURN END SUBROUTINE TDUE(BETA,OMI,THETA,TAUC,NMX,PHI,GAMMA) IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 C(100,100,19) COMMON /A1/ OM(1000,1000),C COMMON /RK10/ SPIN, SI COMMON /A3/ T11,T12,T13 COMMON /TAU1/ TAUS0 COMMON /TAUE/ TAUE COMMON /CONTAT/ ACONT COMMON /MOLFRAZ/ AMOLFRA COMPLEX*16 F0,F1,F2,FMU,FMD COMPLEX*16 A,B,CI,D,EI,F,G,H,AI COMPLEX*16 FF1,FF2,FF3,FF4,FF5,FF6,FF7,FF8,FF9 COMPLEX*16 GG1,GG2,GG3,GG4,GG5,GG6,GG7,GG8,GG9 COMPLEX*16 GEN,FEB,MAR,APR,MAY,JUN,JUL,AUG,SEP COMPLEX*16 QQ1,QQ2,QQ3,QQ4,QQ5,QQ6,QQ7,QQ8,QQ9 COMPLEX*16 TT1,TT2,TT3,TT4,TT5,TT6,TT7,TT8,TT9 COMPLEX*16 HH1,HH2,HH3,HH4,HH5,HH6,HH7,HH8,HH9 COMPLEX*16 AIA,AJ,AT CT=COS(BETA) ST=SIN(BETA) C CONVERT DEG. ---> RAD (CA) CONVER = ATAN(1.0)/45.0 CA=COS(THETA* CONVER) SA=SIN(THETA* CONVER) CF=COS(PHI* CONVER) CF2=COS(2.*PHI* CONVER) SF=SIN(PHI* CONVER) SF2=SIN(2.*PHI* CONVER) C **************** RACAH'S NORMALIZED SPHERICAL HARMONICS F0=-(3.*CA**2-1.)/2. F1=SQRT(6.)/2.*SA*CA*CMPLX(CF,SF) FMU=-SQRT(6.)/2.*SA*CA*CMPLX(CF,-SF) F2=-SQRT(6.)/4.*SA**2*CMPLX(CF2,SF2) FMD=-SQRT(6.)/4.*SA**2*CMPLX(CF2,-SF2) C ELEMENTS OF THE ROTATION MATRIX A=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),SIN(2.*GAMMA)) B=(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA)) CI=(1.+CT)**2/4. D=-(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),SIN(GAMMA)) EI=-ST**2/2. F=-(1.+CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA)) G=(1.-CT)**2/4. H=(1.-CT)*ST/(2.*SQRT(2.))*CMPLX(COS(GAMMA),-SIN(GAMMA)) AI=(1.-CT**2)/4.*CMPLX(COS(2.*GAMMA),-SIN(2.*GAMMA)) GEN=ST**2/2.*CMPLX(COS(2.*GAMMA),SIN(2.*GAMMA)) FEB=CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),SIN(GAMMA)) MAR=-ST**2/2. APR=CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),SIN(GAMMA)) MAY=CT**2 JUN=-CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),-SIN(GAMMA)) JUL=-ST**2/2. AUG=-CT*ST/SQRT(2.)*CMPLX(COS(GAMMA),-SIN(GAMMA)) SEP=ST**2/2.*CMPLX(COS(2.*GAMMA),-SIN(2.*GAMMA)) C ************************ GA=1/TAUC AKL=0. DO 10 K=1,NMX DO 10 L=1,NMX IF(K.EQ.L) GO TO 10 c S+S+ FF1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.) & + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*(3./20.) & + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.)) & *C(K,L,1) c S-S- FF2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.) & + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*(3./20.) & + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.)) & *C(K,L,2) c S-S+ FF3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.) & /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.) & -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(K,L,3) c S+S- FF4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.) & /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.) & -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(K,L,4) c S+SZ FF5=(A*F0*F1*(1./10.)*SQRT(1.5)+B*F0*F0*SQRT(1./50.) & + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+ & G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(K,L,5) c SZS+ FF6=(A*F0*F1*(1./10.)*SQRT(1.5)+D*F0*F0*SQRT(1./50.) & + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+ & CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(K,L,6) c S-SZ FF7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.) & + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+ & CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(K,L,7) c SZS- FF8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.) & + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+ & G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(K,L,8) c SZSZ FF9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI* & F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT(6./50.) & +(3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU) & *C(K,L,9) c S+S+ TT1=(GEN*F0*F0/20. + (FEB+APR)*F0*FMU*(1./20.)*SQRT(3.) & + (MAR+JUL)*F0*FMD*(1./10.)*SQRT(1.5) & + SEP*FMD*FMD*(3./10.) & + MAY*FMU*FMU*(3./20.) & + (JUN+AUG)*FMU*FMD*(3./10.)/SQRT(2.)) & *C(K,L,1) c S-S- TT2=(SEP*F0*F0/20. + (JUN+AUG)*F0*F1*(1./20.)*SQRT(3.) & + (MAR+JUL)*F0*F2*(1./10.)*SQRT(1.5) + MAY*F1*F1*(3./20.) & + (FEB+APR)*F1*F2*(3./10.)/SQRT(2.) + GEN*F2*F2*(3./10.)) & *C(K,L,2) c S-S+ TT3=(-GEN*F0*F2*(1./10.)*SQRT(1.5) & - FEB*F2*FMU*(3./10.)/SQRT(2.) & -MAR*F2*FMD*(3./10.) - APR*F1*F0*(1./20.)*SQRT(3.) & - MAY*F1*FMU*(3./20.) - JUN*F1*FMD*(3./10.)/SQRT(2.) & -JUL*F0*F0*(1./20.) - AUG*F0*FMU*(1./20.)*SQRT(3.) & -SEP*F0*FMD*(1./10.)*SQRT(1.5)) & *C(K,L,3) c S+S- TT4=(-GEN*F0*F2*(1./10.)*SQRT(1.5) & - APR*F2*FMU*(3./10.)/SQRT(2.) & -JUL*F2*FMD*(3./10.) - FEB*F1*F0*(1./20.)*SQRT(3.) & - MAY*F1*FMU*(3./20.) - AUG*F1*FMD*(3./10.)/SQRT(2.) & -MAR*F0*F0*(1./20.) - JUN*F0*FMU*(1./20.)*SQRT(3.) & -SEP*F0*FMD*(1./10.)*SQRT(1.5)) & *C(K,L,4) c S+SZ TT5=(GEN*F0*F1*(1./10.)*SQRT(1.5)+FEB*F0*F0*SQRT(1./50.) & + MAR*F0*FMU*(1./10.)*SQRT(3./2.)+APR*F1*FMU*(3./10.)/SQRT(2.) & +MAY*F0*FMU*SQRT(6./100.)+JUN*FMU*FMU*(3./10.)/SQRT(2.)+ & JUL*F1*FMD*3./10.+AUG*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.) & *C(K,L,5) c SZS+ TT6=(GEN*F0*F1*(1./10.)*SQRT(1.5)+APR*F0*F0*SQRT(1./50.) & + JUL*F0*FMU*(1./10.)*SQRT(3./2.)+FEB*F1*FMU*(3./10.)/SQRT(2.) & +MAY*F0*FMU*SQRT(6./100.)+AUG*FMU*FMU*(3./10.)/SQRT(2.)+ & MAR*F1*FMD*3./10.+JUN*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.) & *C(K,L,6) c S-SZ TT7=-(SEP*F0*FMU*(1./10.)*SQRT(1.5)+AUG*F0*F0*SQRT(1. & /50.)+JUL*F0*F1*(1./10.)*SQRT(3./2.)+JUN*F1*FMU*(3./10.)/SQRT(2.) & +MAY*F0*F1*SQRT(6./100.)+APR*F1*F1*(3./10.)/SQRT(2.)+ & MAR*FMU*F2*3./10.+FEB*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.) & *C(K,L,7) c SZS- TT8=-(SEP*F0*FMU*(1./10.)*SQRT(1.5)+JUN*F0*F0* & SQRT(1./50.)+MAR*F0*F1*(1./10.)*SQRT(3./2.)+AUG*F1*FMU*(3./10.) & /SQRT(2.)+MAY*F0*F1*SQRT(6./100.)+FEB*F1*F1*(3./10.)/SQRT(2.)+ & JUL*FMU*F2*3./10.+APR*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.) & *C(K,L,8) c SZSZ TT9=((3./10.)*GEN*F1*F1+FEB*F1*F0*SQRT(6./50.) & +(3./10.)*MAR*F1*FMU+ & APR*F0*F1*SQRT(6./50)+(2./5.)*MAY*F0*F0+ & JUN*F0*FMU*SQRT(6./50.)+ & (3./10.)*JUL*F1*FMU+AUG*FMU*F0*SQRT(6./50.) & +(3./10.)*SEP*FMU*FMU) & *C(K,L,9) c S+S+ GG1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.) & + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*(3./20.) & + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.)) & *C(L,K,1) c S-S- GG2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.) & + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*(3./20.) & + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.)) & *C(L,K,2) c S-S+ GG3=(-A*F0*F2*(1./10.)*SQRT(1.5) - B*F2*FMU*(3./10.) & /SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.) & -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(L,K,3) c S+S- GG4=(-A*F0*F2*(1./10.)*SQRT(1.5) - D*F2*FMU*(3./10.) & /SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.) & -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(L,K,4) c S+SZ GG5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.) & + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+ & G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(L,K,5) c SZS+ GG6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.) & + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+ & CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(L,K,6) c S-SZ GG7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.) & + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+ & CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(L,K,7) c SZS- GG8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.) & + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+ & G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(L,K,8) c SZSZ GG9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)* & CI*F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT(6./ & 50.)+(3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU) & *C(L,K,9) AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX) AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX) AT=(TT1+TT2+TT3+TT4+TT5+TT6+TT7+TT8+TT9)/FLOAT(NMX) TMP=0.5*(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))* & 1.E+9/((OM(K,L)+OMI)**2+GA**2)+ & 0.5*(REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))* & 1.E+9/((OM(L,K)-OMI)**2+GA**2) & -(REAL(AT)*GA+IMAG(AT)*OM(K,L))* & 1.E+9/(OM(K,L)**2+GA**2) AKL=TMP+AKL 10 CONTINUE c S+S+ HH1=(A*F0*F0/20. + (B+D)*F0*FMU*(1./20.)*SQRT(3.) & + (CI+G)*F0*FMD*(1./10.)*SQRT(1.5) + EI*FMU*FMU*(3./20.) & + (F+H)*FMU*FMD*(3./10.)/SQRT(2.) + AI*FMD*FMD*(3./10.)) & *C(1,1,11) c S-S- HH2=(AI*F0*F0/20. + (F+H)*F0*F1*(1./20.)*SQRT(3.) & + (CI+G)*F0*F2*(1./10.)*SQRT(1.5) + EI*F1*F1*(3./20.) & + (B+D)*F1*F2*(3./10.)/SQRT(2.) + A*F2*F2*(3./10.)) & *C(1,1,12) c S-S+ HH3=(-A*F0*F2*(1./10.)*SQRT(1.5)-B*F2*FMU*(3./10.)/ & SQRT(2.)-CI*F2*FMD*(3./10.) - D*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - F*F1*FMD*(3./10.)/SQRT(2.) & -G*F0*F0*(1./20.) - H*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(1,1,13) c S+S- HH4=(-A*F0*F2*(1./10.)*SQRT(1.5)-D*F2*FMU*(3./10.)/ & SQRT(2.)-G*F2*FMD*(3./10.) - B*F1*F0*(1./20.)*SQRT(3.) & - EI*F1*FMU*(3./20.) - H*F1*FMD*(3./10.)/SQRT(2.) & -CI*F0*F0*(1./20.) - F*F0*FMU*(1./20.)*SQRT(3.) & -AI*F0*FMD*(1./10.)*SQRT(1.5)) & *C(1,1,14) c S+SZ HH5=(A*F0*F1*(1./10.)*SQRT(1.5) + B*F0*F0*SQRT(1./50.) & + CI*F0*FMU*(1./10.)*SQRT(3./2.)+D*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+F*FMU*FMU*(3./10.)/SQRT(2.)+ & G*F1*FMD*3./10.+H*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(1,1,15) c SZS+ HH6=(A*F0*F1*(1./10.)*SQRT(1.5) + D*F0*F0*SQRT(1./50.) & + G*F0*FMU*(1./10.)*SQRT(3./2.)+B*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*FMU*SQRT(6./100.)+H*FMU*FMU*(3./10.)/SQRT(2.)+ & CI*F1*FMD*3./10.+F*F0*FMD*SQRT(6./50.)+AI*FMU*FMD*3./10.) & *C(1,1,16) c S-SZ HH7=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+H*F0*F0*SQRT(1./50.) & + G*F0*F1*(1./10.)*SQRT(3./2.)+F*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+D*F1*F1*(3./10.)/SQRT(2.)+ & CI*FMU*F2*3./10.+B*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(1,1,17) c SZS- HH8=-(AI*F0*FMU*(1./10.)*SQRT(1.5)+F*F0*F0*SQRT(1./50.) & + CI*F0*F1*(1./10.)*SQRT(3./2.)+H*F1*FMU*(3./10.)/SQRT(2.)+ & EI*F0*F1*SQRT(6./100.)+B*F1*F1*(3./10.)/SQRT(2.)+ & G*FMU*F2*3./10.+D*F0*F2*SQRT(6./50.)+A*F1*F2*3./10.) & *C(1,1,18) c SZSZ HH9=((3./10.)*A*F1*F1+B*F1*F0*SQRT(6./50.)+(3./10.)*CI* & F1*FMU+D*F0*F1*SQRT(6./50)+(2./5.)*EI*F0*F0+F*F0*FMU*SQRT(6./50. & )+(3./10.)*G*F1*FMU+H*FMU*F0*SQRT(6./50.)+(3./10.)*AI*FMU*FMU) & *C(1,1,19) c S+S+ QQ1=(GEN*F0*F0/20. + (FEB+APR)*F0*FMU*(1./20.)*SQRT(3.) & + (MAR+JUL)*F0*FMD*(1./10.)*SQRT(1.5) & + SEP*FMD*FMD*(3./10.) & + MAY*FMU*FMU*(3./20.) & + (JUN+AUG)*FMU*FMD*(3./10.)/SQRT(2.)) & *C(1,1,11) c S-S- QQ2=(SEP*F0*F0/20. + (JUN+AUG)*F0*F1*(1./20.)*SQRT(3.) & + (MAR+JUL)*F0*F2*(1./10.)*SQRT(1.5) + MAY*F1*F1*(3./20.) & + (FEB+APR)*F1*F2*(3./10.)/SQRT(2.) + GEN*F2*F2*(3./10.)) & *C(1,1,12) c S-S+ QQ3=(-GEN*F0*F2*(1./10.)*SQRT(1.5) & - FEB*F2*FMU*(3./10.)/SQRT(2.) & -MAR*F2*FMD*(3./10.) - APR*F1*F0*(1./20.)*SQRT(3.) & - MAY*F1*FMU*(3./20.) - JUN*F1*FMD*(3./10.)/SQRT(2.) & -JUL*F0*F0*(1./20.) - AUG*F0*FMU*(1./20.)*SQRT(3.) & -SEP*F0*FMD*(1./10.)*SQRT(1.5)) & *C(1,1,13) c S+S- QQ4=(-GEN*F0*F2*(1./10.)*SQRT(1.5) & - APR*F2*FMU*(3./10.)/SQRT(2.) & -JUL*F2*FMD*(3./10.) - FEB*F1*F0*(1./20.)*SQRT(3.) & - MAY*F1*FMU*(3./20.) - AUG*F1*FMD*(3./10.)/SQRT(2.) & -MAR*F0*F0*(1./20.) - JUN*F0*FMU*(1./20.)*SQRT(3.) & -SEP*F0*FMD*(1./10.)*SQRT(1.5)) & *C(1,1,14) c S+SZ QQ5=(GEN*F0*F1*(1./10.)*SQRT(1.5) + FEB*F0*F0*SQRT( & 1./50.)+MAR*F0*FMU*(1./10.)*SQRT(3./2.)+APR*F1*FMU*(3./10.)/ & SQRT(2.)+MAY*F0*FMU*SQRT(6./100.)+JUN*FMU*FMU*(3./10.)/SQRT(2.)+ & JUL*F1*FMD*3./10.+AUG*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.) & *C(1,1,15) c SZS+ QQ6=(GEN*F0*F1*(1./10.)*SQRT(1.5) + APR*F0*F0*SQRT & (1./50.)+JUL*F0*FMU*(1./10.)*SQRT(3./2.)+FEB*F1*FMU*(3./10.)/ & SQRT(2.)+MAY*F0*FMU*SQRT(6./100.)+AUG*FMU*FMU*(3./10.)/SQRT(2.)+ & MAR*F1*FMD*3./10.+JUN*F0*FMD*SQRT(6./50.)+SEP*FMU*FMD*3./10.) & *C(1,1,16) c S-SZ QQ7=-(SEP*F0*FMU*(1./10.)*SQRT(1.5) + AUG*F0*F0*SQRT & (1./50.)+ JUL*F0*F1*(1./10.)*SQRT(3./2.)+JUN*F1*FMU*(3./10.)/ & SQRT(2.)+MAY*F0*F1*SQRT(6./100.)+APR*F1*F1*(3./10.)/SQRT(2.)+ & MAR*FMU*F2*3./10.+FEB*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.) & *C(1,1,17) c SZS- QQ8=-(SEP*F0*FMU*(1./10.)*SQRT(1.5) + JUN*F0*F0*SQRT & (1./50.)+ MAR*F0*F1*(1./10.)*SQRT(3./2.)+AUG*F1*FMU*(3./10.)/ & SQRT(2.)+MAY*F0*F1*SQRT(6./100.)+FEB*F1*F1*(3./10.)/SQRT(2.)+ & JUL*FMU*F2*3./10.+APR*F0*F2*SQRT(6./50.)+GEN*F1*F2*3./10.) & *C(1,1,18) c SZSZ QQ9=((3./10.)*GEN*F1*F1+FEB*F1*F0*SQRT(6./50.) & +(3./10.)*MAR*F1*FMU+ & APR*F0*F1*SQRT(6./50)+(2./5.)*MAY*F0*F0+ & JUN*F0*FMU*SQRT(6./50.)+ & (3./10.)*JUL*F1*FMU+AUG*FMU*F0*SQRT(6./50.) & +(3./10.)*SEP*FMU*FMU) & *C(1,1,19) ALL=(REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX)) & *GA*1.E+9/(OMI**2+GA**2)+(IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7 & +HH8+HH9)/FLOAT(NMX)*OMI-IMAG(HH1+HH2+HH3+HH4+HH5+HH6+HH7+ & HH8+HH9)/FLOAT(NMX)*OMI) & *1.E+9/(OMI**2+GA**2) & -(REAL(QQ1+QQ2+QQ3+QQ4+QQ5+QQ6+QQ7+QQ8+QQ9)/FLOAT(NMX))* & 1.E+9*(1./GA) C DIPOLAR TERM T11=(-ALL-AKL) C CONTACT TERM IF (ACONT.NE.0)THEN IF(AMOLFRA.EQ.0)THEN RKCONT=1./(SPIN*(SPIN+1.)*2./3./ & (1.0546)**2*1.E34*1.E34* & (ACONT*6.28*1.0546E-34*1.E6)**2)/TAUS0 ELSE RKCONT=1. ENDIF C **************** F0=(1.) F1=0. F2=0. FMU=0. FMD=0. C ************************ GA=1/TAUE T1CONT=0. DO 101 K=1,NMX DO 101 L=1,NMX IF(K.EQ.L) GO TO 101 c S+S+ FF1=(A*F0*F0/2.) & *C(K,L,1) c S-S- FF2=(AI*F0*F0/2) & *C(K,L,2) c S+S- FF3=-G*F0*F0/2. & *C(K,L,3) c S-S+ FF4=-CI*F0*F0/2 & *C(K,L,4) c SZS+ FF5=-B*F0*F0*SQRT(1./2.) & *C(K,L,5) c S+SZ FF6=-D*F0*F0*SQRT(1./2.) & *C(K,L,6) c SZS- FF7=H*F0*F0*SQRT(1./2.) & *C(K,L,7) c S-SZ FF8=F*F0*F0*SQRT(1./2.) & *C(K,L,8) c SZSZ FF9=EI*F0*F0 & *C(K,L,9) c S+S+ GG1=(A*F0*F0/2.) & *C(L,K,1) c S-S- GG2=(AI*F0*F0/2) & *C(L,K,2) c S+S- GG3=-G*F0*F0/2. & *C(L,K,3) c S-S+ GG4=-CI*F0*F0/2 & *C(L,K,4) c SZS+ GG5=-B*F0*F0*SQRT(1./2.) & *C(L,K,5) c S+SZ GG6=-D*F0*F0*SQRT(1./2.) & *C(L,K,6) c SZS- GG7=H*F0*F0*SQRT(1./2.) & *C(L,K,7) c S-SZ GG8=F*F0*F0*SQRT(1./2.) & *C(L,K,8) c SZSZ GG9=EI*F0*F0 & *C(L,K,9) c S+S+ TT1=GEN/2.*F0*F0*C(K,L,1) c S-S- TT2=SEP/2.*F0*F0*C(K,L,2) c S-S+ TT3=-JUL/2.*F0*F0 & *C(K,L,3) c S+S- TT4=-MAR*F0*F0/2 & *C(K,L,4) c S+SZ TT5=-FEB*F0*F0/SQRT(2.) & *C(K,L,5) c SZS+ TT6=-APR*F0*F0/SQRT(2.) & *C(K,L,6) c S-SZ TT7=AUG*F0*F0/SQRT(2.) & *C(K,L,7) c SZS- TT8=JUN*F0*F0/SQRT(2.)*C(K,L,8) c SZSZ TT9=MAY*F0*F0 & *C(K,L,9) AIA=(FF1+FF2+FF3+FF4+FF5+FF6+FF7+FF8+FF9)/FLOAT(NMX) AJ=(GG1+GG2+GG3+GG4+GG5+GG6+GG7+GG8+GG9)/FLOAT(NMX) AT=(TT1+TT2+TT3+TT4+TT5+TT6+TT7+TT8+TT9)/FLOAT(NMX) TMP=0.5*(REAL(AIA)*GA+IMAG(AIA)*(OMI+OM(K,L)))* & 1.E+34/((OM(K,L)+OMI)**2+GA**2)+ & 0.5*(REAL(AJ)*GA+IMAG(AJ)*(-OMI+OM(L,K)))* & 1.E+34/((OM(L,K)-OMI)**2+GA**2) & -(REAL(AT)+IMAG(AT)*OM(K,L))*1.E+34*GA/(OM(K,L)**2+GA**2) T1CONT=-TMP+T1CONT 101 CONTINUE c S+S+ HH1=(A*F0*F0/2.) & *C(1,1,11) c S-S- HH2=(AI*F0*F0/2) & *C(1,1,12) c S+S- HH3=-G*F0*F0/2. & *C(1,1,13) c S-S+ HH4=-CI*F0*F0/2 & *C(1,1,14) c SZS+ HH5=-B*F0*F0*SQRT(1./2.) & *C(1,1,15) c S+SZ HH6=-D*F0*F0*SQRT(1./2.) & *C(1,1,16) c SZS- HH7=H*F0*F0*SQRT(1./2.) & *C(1,1,17) c S-SZ HH8=F*F0*F0*SQRT(1./2.) & *C(1,1,18) c SZSZ HH9=EI*F0*F0 & *C(1,1,19) c S+S+ QQ1=GEN/2.*F0*F0*C(1,1,11) c S-S- QQ2=SEP/2.*F0*F0*C(1,1,12) c S-S+ QQ3=-JUL/2.*F0*F0 & *C(1,1,13) c S+S- QQ4=-MAR*F0*F0/2. & *C(1,1,14) c S+SZ QQ5=-FEB*F0*F0/SQRT(2.) & *C(1,1,15) c SZS+ QQ6=-APR*F0*F0/SQRT(2.) & *C(1,1,16) c S-SZ QQ7=AUG*F0*F0/SQRT(2.) & *C(1,1,17) c SZS- QQ8=JUN*F0*F0/SQRT(2.)*C(1,1,18) c SZSZ QQ9=MAY*F0*F0 & *C(1,1,19) C TERMS FROM INT(EXP(-T/TAU)DT),SZ-LABORATORY COORD.FRAME TMQ=REAL(QQ1+QQ2+QQ3+QQ4+QQ5+QQ6+QQ7+QQ8+QQ9)/FLOAT(N) & *1.E+34*(1./GA) C TERMS FROM INT(EXP(IOMI*T-T*GA)DT + C EXP(-OMI*T-T*GA)DT),S+,S- - LABORATORY COORD.FRAME TH=2*REAL(HH1+HH2+HH3+HH4+HH5+HH6+HH7+HH8+HH9)/FLOAT(NMX)*GA & *1.E+34*GA/(OMI**2+GA**2) C CONTACT CONTRIBUTION T1CONT=(T1CONT-TH+TMQ)/(1.0546)**2*1.E34 T12=T1CONT*RKCONT ENDIF RETURN END SUBROUTINE TUNOISO(BETA,OMI,THETA,TAUC,NMX) IMPLICIT REAL*8(A-H,O-Z) C CALCOLA IVALORI DI T1**-1 CHE PERO* DEVONO ESSERE ANCORA INTEGRATI S COMMON /AOLD/ OM(10000),C(10000,4) COMMON /TAUE/ TAUE COMMON /CONTAT/ ACONT COMMON /TAU1/ TAUS0 COMMON /RK10/ SPIN, SI COMMON /MOLFRAZ/ AMOLFRA COMMON /A3/ T11,T12,T13 JX(I)=(NMX-1)*I-(NMX-3)-NMX*(NMX-2)*((IP-2)/NMX) CT=DCOS(BETA) ST=DSIN(BETA) C CONVERT DEG. ---> RAD (CA) CONVER = ATAN(1.0)/45.0 C **************** CA=DCOS(THETA* CONVER) FZ=(3.*CA**2-1.)**2 FU=9.*CA**2*(1.-CA**2)/4. FD=9.*(1.-CA**2)**2/16. C ELEMENTS OF THE ROTATION MATRIX A=-(1.-CT**2)/4. B=-(1.+CT)*ST/4. CI=(1.+CT)**2/4. D=(1.-CT)*ST/4. EI=ST**2/4. F=(1.-CT)**2/4. C ************************ GA=1/TAUC AKL=0. IP=2 IA=NMX 20 CONTINUE DO 10 I=IP,IA J=JX(I) F1=A*FZ/8.*(C(J,1)+C(I,1)) F2=-(B*FZ-4.*D*FU)*C(I,2) F3=-(D*FZ-4.*B*FU)*C(J,2) F4=(CI*FZ/8.+2.*EI*FU+2.*F*FD)*C(I,3) F5=(F*FZ/8.+2.*EI*FU+2.*CI*FD)*C(J,3) F6=2.*(EI*FZ+CI*FU)*C(I,4)+2.*F*FU*C(J,4) G1=A*FZ/8.*(C(J,1)+C(I,1)) G2=-(B*FZ-4.*D*FU)*C(J,2) G3=-(D*FZ-4.*B*FU)*C(I,2) G4=(CI*FZ/8.+2.*EI*FU+2.*F*FD)*C(J,3) G5=(F*FZ/8.+2.*EI*FU+2.*CI*FD)*C(I,3) G6=2.*(EI*FZ+CI*FU)*C(J,4)+2.*F*FU*C(I,4) AI=(F1+F2+F3+F4+F5+F6)/FLOAT(NMX) AJ=(G1+G2+G3+G4+G5+G6)/FLOAT(NMX) TMP=2.*(AI+AJ)*GA*1.E+9/(OM(I)**2+GA**2) AKL=TMP+AKL 10 CONTINUE IP=IP+NMX IA=IA+NMX-1 IF (IA.LT.IP) GO TO 30 GO TO 20 30 CONTINUE H1=2.*A*FZ/8.*C(1,1) H2=-(B*FZ-4.*D*FU)*C(1,2) H3=-(D*FZ-4.*B*FU)*C(1,2) H4=(CI*FZ/8.+2.*EI*FU+2.*F*FD)*C(1,3) H5=(F*FZ/8.+2.*EI*FU+2.*CI*FD)*C(1,3) H6=2.*(EI*FZ+CI*FU)*C(1,4)+2.*F*FU*C(1,4) ALL=2.*(H1+H2+H3+H4+H5+H6)/FLOAT(NMX)*GA*1.E+9/(OMI**2+GA**2) C DIPOLAR TERM T11=(AKL+ALL)/10. IF (ACONT.NE.0)THEN IF(AMOLFRA.EQ.0)THEN RKCONT=1./(SPIN*(SPIN+1.)*2./3./ & (1.0546)**2*1.E34*1.E34* & (ACONT*6.28*1.0546E-34*1.E6)**2)/TAUS0 ELSE RKCONT=1. ENDIF CT=DCOS(BETA) ST=DSIN(BETA) C CONVERT DEG. ---> RAD (CA) CONVER = ATAN(1.0)/45.0 C **************** CA=DCOS(THETA* CONVER) FZ=(1.) FU=0. FD=0. C ELEMENTS OF THE ROTATION MATRIX A=-(1.-CT**2)/4. B=-(1.+CT)*ST/4. CI=(1.+CT)**2/4. D=(1.-CT)*ST/4. EI=ST**2/4. F=(1.-CT)**2/4. C ************************ GA=1/TAUE T1CONT=0. IP=2 IA=NMX 120 CONTINUE DO 110 I=IP,IA J=JX(I) F1=A*FZ/2.*(C(J,1)+C(I,1)) F2=-(-B*FZ*2-4.*D*FU)*C(I,2) F3=-(-D*FZ*2-4.*B*FU)*C(J,2) F4=(CI*FZ/2.+2.*EI*FU+2.*F*FD)*C(I,3) F5=(F*FZ/2.+2.*EI*FU+2.*CI*FD)*C(J,3) F6=2.*(EI*FZ+CI*FU)*C(I,4)+2.*F*FU*C(J,4) G1=A*FZ/2.*(C(J,1)+C(I,1)) G2=-(-B*FZ*2-4.*D*FU)*C(J,2) G3=-(-D*FZ*2-4.*B*FU)*C(I,2) G4=(CI*FZ/2.+2.*EI*FU+2.*F*FD)*C(J,3) G5=(F*FZ/2.+2.*EI*FU+2.*CI*FD)*C(I,3) G6=2.*(EI*FZ+CI*FU)*C(J,4)+2.*F*FU*C(I,4) AI=(F1+F2+F3+F4+F5+F6)/FLOAT(NMX) AJ=(G1+G2+G3+G4+G5+G6)/FLOAT(NMX) TMP=2.*(AI+AJ)*GA*1.E34/(OM(I)**2+GA**2) C CONTACT TERM T1CONT=TMP+T1CONT 110 CONTINUE IP=IP+NMX IA=IA+NMX-1 IF (IA.LT.IP) GO TO 130 GO TO 120 130 CONTINUE T12=T1CONT/(1.0546)**2*1.E34*RKCONT ENDIF RETURN END SUBROUTINE suPOWELL(P,XI,N,NP,FTOL,ITER,FRET,NMX) C PERFORMES THE FITTING PROCEDURE IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=20,ITMAX=2000) DIMENSION P(NP),XI(NP,NP),PT(NMAX),PTT(NMAX),XIT(NMAX) COMMON /NMXCOM/ NMX1 COMMON /NVNP/ NV NMX1=NMX NV=NP CALL XISTEP(NV,XI,P) CALL FUNCsuper(P,FUNC,NMX,NV) FRET=FUNC DO 11 J=1,N PT(J)=P(J) 11 CONTINUE ITER=0 1 ITER=ITER+1 FP=FRET IBIG=0 DEL=0. DO 13 I=1,N DO 12 J=1,N XIT(J)=XI(J,I) 12 CONTINUE FPTT=FRET CALL LINMIN3(P,XIT,N,FRET) IF(ABS(FPTT-FRET).GT.DEL)THEN DEL=ABS(FPTT-FRET) IBIG=I ENDIF 13 CONTINUE IF(2.*ABS(FP-FRET).LE.FTOL*(ABS(FP)+ABS(FRET)))RETURN IF(ITER.EQ.ITMAX) PAUSE 'POWELL EXCEEDING MAXIMUM ITERATIONS.' DO 14 J=1,N PTT(J)=2.*P(J)-PT(J) XIT(J)=P(J)-PT(J) PT(J)=P(J) 14 CONTINUE CALL FUNCsuper(P,FUNC,NMX,NV) FPTT=FUNC IF(FPTT.GE.FP)GO TO 1 T=2.*(FP-2.*FRET+FPTT)*(FP-FRET-DEL)**2-DEL*(FP-FPTT)**2 IF(T.GE.0.)GO TO 1 CALL LINMIN3(P,XIT,N,FRET) DO 15 J=1,N XI(J,IBIG)=XIT(J) 15 CONTINUE GO TO 1 END FUNCTION F3DIM(X) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=50) COMMON /F3COM/ PCOM(NMAX),XICOM(NMAX),NCOM COMMON /NMXCOM/ NMX1 COMMON /NVNP/ NV DIMENSION XT(NMAX) DO 11 J=1,NCOM XT(J)=PCOM(J)+X*XICOM(J) 11 CONTINUE CALL FUNCsuper(xt,FUNC,NMX1,NV) F3DIM=FUNC RETURN END SUBROUTINE LINMIN3(P,XI,N,FRET) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=50,TOL=1.E-4) EXTERNAL F3DIM DIMENSION P(N),XI(N) COMMON /F3COM/ PCOM(NMAX),XICOM(NMAX),NCOM NCOM=N DO 11 J=1,N PCOM(J)=P(J) XICOM(J)=XI(J) 11 CONTINUE AX=0. XX=1 ! BX=2 CALL MNBRAK3(AX,XX,BX,FA,FX,FB,F3DIM) FRET=BRENT3(AX,XX,BX,F3DIM,TOL,XMIN) DO 12 J=1,N XI(J)=XMIN*XI(J) P(J)=P(J)+XI(J) 12 CONTINUE RETURN END SUBROUTINE MNBRAK3(AX,BX,CX,FA,FB,FC,F3DIM) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20) FA=F3DIM(AX) FB=F3DIM(BX) IF(FB.GT.FA)THEN DUM=AX AX=BX BX=DUM DUM=FB FB=FA FA=DUM ENDIF CX=BX+GOLD*(BX-AX) FC=F3DIM(CX) 1 IF(FB.GE.FC)THEN R=(BX-AX)*(FB-FC) Q=(BX-CX)*(FB-FA) U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R)) ULIM=BX+GLIMIT*(CX-BX) IF((BX-U)*(U-CX).GT.0.)THEN FU=F3DIM(U) IF(FU.LT.FC)THEN AX=BX FA=FB BX=U FB=FU GO TO 1 ELSE IF(FU.GT.FB)THEN CX=U FC=FU GO TO 1 ENDIF U=CX+GOLD*(CX-BX) FU=F3DIM(U) ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN FU=F3DIM(U) IF(FU.LT.FC)THEN BX=CX CX=U U=CX+GOLD*(CX-BX) FB=FC FC=FU FU=F3DIM(U) ENDIF ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN U=ULIM FU=F3DIM(U) ELSE U=CX+GOLD*(CX-BX) FU=F3DIM(U) ENDIF AX=BX BX=CX CX=U FA=FB FB=FC FC=FU GO TO 1 ENDIF RETURN END FUNCTION BRENT3(AX,BX,CX,F3DIM,TOL,XMIN) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10) A=MIN(AX,CX) B=MAX(AX,CX) V=BX W=V X=V E=0. FX=F3DIM(X) FV=FX FW=FX DO 11 ITER=1,ITMAX XM=0.5*(A+B) TOL1=TOL*ABS(X)+ZEPS TOL2=2.*TOL1 IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3 IF(ABS(E).GT.TOL1) THEN R=(X-W)*(FX-FV) Q=(X-V)*(FX-FW) P=(X-V)*Q-(X-W)*R Q=2.*(Q-R) IF(Q.GT.0.) P=-P Q=ABS(Q) ETEMP=E E=D IF(ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR. * P.GE.Q*(B-X)) GOTO 1 D=P/Q U=X+D IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X) GOTO 2 ENDIF 1 IF(X.GE.XM) THEN E=A-X ELSE E=B-X ENDIF D=CGOLD*E 2 IF(ABS(D).GE.TOL1) THEN U=X+D ELSE U=X+SIGN(TOL1,D) ENDIF FU=F3DIM(U) IF(FU.LE.FX) THEN IF(U.GE.X) THEN A=X ELSE B=X ENDIF V=W FV=FW W=X FW=FX X=U FX=FU ELSE IF(U.LT.X) THEN A=U ELSE B=U ENDIF IF(FU.LE.FW .OR. W.EQ.X) THEN V=W FV=FW W=U FW=FU ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN V=U FV=FU ENDIF ENDIF 11 CONTINUE ! PAUSE 'BRENT EXCEED MAXIMUM ITERATIONS.' WRITE(6,*) 'BRENT3 EXCEED MAXIMUM ITERATIONS.' RETURN !!!!!!!!!!!!!!!!!!!!! 3 XMIN=X BRENT3=FX RETURN END SUBROUTINE provag(chbo,out) common /super/ b,DIF,xMSAT,TAUN,HANI,TEMP INTEGER ji INTEGER j,j2 REAL B,DIF COMPLEX I,R integer mbranc, nmx PARAMETER(MBRANC=21) ! DIMENSION WR( MBRANC ) ! DIMENSION ARR(MBRANC,MBRANC),ARI(MBRANC,MBRANC) ! DIMENSION ZRR(MBRANC,MBRANC),ZRI(MBRANC,MBRANC) ! DIMENSION WK1(MBRANC),WK2(MBRANC),WK3(MBRANC) real*8 WR( MBRANC ) real*8 ARR(MBRANC,MBRANC),ARI(MBRANC,MBRANC) real*8 ZRR(MBRANC,MBRANC),ZRI(MBRANC,MBRANC) real*8 WK1(MBRANC),WK2(MBRANC),WK3(MBRANC) REAL, DIMENSION (90) :: mdensangl REAL, DIMENSION (101) :: dd,tt,ddd REAL, DIMENSION (2,21) :: A REAL, DIMENSION (21) :: eval,mdens,r1,r2 REAL, DIMENSION (21,21) :: evec,invec,tempon REAL, DIMENSION (21,21):: SAPL,SANEG,SAZ, & SZZD,SZPLD,SZNEGD,jnegze,jzze REAL, DIMENSION (21,21) :: SZZ ,JZAYANT,SZPL,SZNEG OPEN(UNIT=3,FILE='spin50.dat') OPEN(UNIT=4,FILE='spin50.') ! hskghz:rapport entre c Plank et c Boltzman si frequence GHz ! do 7000 nex=1,4 ! hani=ddd(nex) ! b=dd(nex) Write(4,*) "courbe",nex,"hani:",hani,"d:",b Write(3,*)"courbe",nex,"hani:",hani,"d:",b ! Write(*,*) "courbe",nex,"hani:",hani,"d:",b hskghz=0.06626/1.3181 ! spinmul ratio between real and fictif spin TCD=(B*B)/dif xMASPAR=3.1416*4*B*B*B*5.18/3 sisatu=xmsat*xmaspar/1000 nelec=sisatu/9.274e-24 spinmul=nelec/2 ! write(*,*) nspin,spinmul xMSATU=xMSAT*xMASPAR CONCPAR=232E-3/(3*xMASPAR) CMT=1.7765e+5*CONCPAR*xMSATU*xMSATU/(b*dif) Spin=10 xspin=10. Spin2=xSpin*(xSpin+1.) CMT=cmt/spin2 Nspin=INT(2*Spin+1) r1=0. r2=0. SAPL=0. SANEG=0. SAZ=0. ! Calculation in the anisotropy axis system Ba DO ni=1,20 saz(ni,ni)=real(spin-ni+1) sapl(ni,ni+1)=sqrt(Real(spin2-(spin-ni+1)*(spin-ni))) saneg(ni+1,ni)=sqrt(real(spin2-(spin-ni+1)*(spin-ni))) end do saz(21,21)=-spin n=21 lda=2 ncoda=1 ldevec=21 pasbo=1.123 !in the routine the energy is express in GHz; chbo is the external field ! chbo=0.00001/pasbo somnor=0. ! DO 5000 nn=1,100,1 ! CHBO=CHBO*pasbo r1(1)=0. r1(2)=0. r1(3)=0. r2(1)=0. r2(2)=0. r2(3)=0. sometang=0.0 do 6000 nagl=0,90,1 angl=1.*nagl !Expression of Sz,S+,S- operator in anisotropy axes system SZZ=cosd(angl)*SAZ+0.5*sind(angl)*(SAPL+SANEG) SZPL=0.5*cosd(angl)*(SAPL+SANEG)+0.5*(SAPL-SANEG)- & sind(angl)*SAZ SZNEG=+0.5*cosd(angl)*(SAPL+SANEG)-0.5*(SAPL-SANEG)- & sind(angl)*SAZ EZEEMAN=CHBO*658. OMEGA=CHBO*1000000000. consa1=hani/(2.*xSpin) consa2=EZEEMAN*cosd(angl) consc=0.5*EZEEMAN*sind(angl) do ni=1,Nspin Smz=xSpin+1.-ni A(2,ni)=-(consa1*Smz*Smz+consa2*Smz) ! write(6,*)'a',a(2,ni) end do A(1,1)=0 do ni=2,Nspin Smz=xSpin+1.-ni A(1,ni)=-consc*SQRT(real(Spin2-Smz*(Smz+1))) ! write(6,*)'a',a(1,ni) end do ! calculation of the eight and vector value C eval EIGENVALUES evec EIGENVECTORS do ni=1,Nspin ARR(ni,ni)=A(2,ni) end do do ni=2,Nspin ARR(ni,ni-1)=A(1,ni) end do do ni=1,Nspin-1 ARR(ni,ni+1)=A(1,ni) end do ! write(6,*)' ' ! DO 247 J2=1,Nspin ! write(6,*)arr(J2,J2) ! DO 247 J=1,Nspin !247 CONTINUE ! write(6,*)' ' DO 245 J2=1,ni DO 245 J=1,ni ARI(J2,J)=0. 245 CONTINUE CALL F02AXF(ARR,MBRANC,ARI,MBRANC,n,wr,zrr,MBRANC, & ZRI,MBRANC,WK1,WK2,WK3,0) ! write(6,*)Nspin DO 246 J2=1,Nspin eval(J2)=wr(J2) DO 246 J=1,Nspin evec(J,J2)=zrr(J,J2) 246 CONTINUE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 ! CALL EVCSB(N,A,LDA,NCODA,EVAL,EVEC,LDEVEC) do ni=1,Nspin do nj=1,Nspin invec(nj,ni)=evec(ni,nj) end do end do ! Calculation of the diagonal term of the density matrix and normalization xnor1=0.0 sumener=0.0 expoa= ABS(-spinmul*hskghz*eval(1)/(temp*spin)) ! Procedure to avoid probem with exponent to large n10=INT(expoa/10)-1 do ns=1,Nspin expob=-spinmul*hskghz*eval(ns)/(temp*spin) expoc=expob-n10*10. mdens(ns)=0. if (expoc.gt.-30.0) then mdens(ns)=exp(expoc) end if xnor1=xnor1+mdens(ns) sumener=sumener+mdens(ns)*eval(ns) end do do nde=1,nspin mdens(nde)=mdens(nde)/xnor1 end do enmoy=sumener/xnor1 if (nagl.eq.0) then exangl= -spinmul*hskghz*enmoy/(temp*spin) ! Procedure to avoid probem with exponent to large end if expanglb=-spinmul*hskghz*enmoy/(temp*spin) expanglc=expanglb-exangl mdensangl(nagl)=0. if (expanglc.gt.-10) then mdensangl(nagl)=sind(angl)*exp(expanglc) sometang=sometang+mdensangl(nagl) end if tempon=MATMUL(szz,evec) SZZD=MATMUL(invec,tempon) tempon=MATMUL(szpl,evec) SZplD=MATMUL(invec,tempon) tempon=MATMUL(szneg,evec) SZnegd=MATMUL(invec,tempon) Szm=0 do ns=1,Nspin Szm=Szzd(ns,ns)*mdens(ns)+Szm end do szm2=szm*szm omegai=1000000000.*chbo omegaii=6.28*omegai rzxfreed = 0 rxxfreed = 0 rzzfreed = 0 rxzfreed = 0 rayan1=0. rayan2=0. raya=rayant(cmt,omegaii,tcd) rayz=rayant(cmt,0.,tcd) do ni=1,Nspin taun1=taun*exp(spinmul*hskghz*(eval(1)-eval(ni))/(temp*spin)) if(taun1.lt.1e-16)taun1=1e-16 do nj=1,Nspin omegae=1000000000.*(eval(ni)-eval(nj)) omega=6.28*abs(omegae+omegai) omegan=6.28*abs(omegae-omegai) omegaz=6.28*abs(omegae) rfie=rfreed(cmt,omega,taun1,tcd) rfien=rfreed(cmt,omegan,taun1,tcd) rfe=rfreed(cmt,omegaz,taun1,tcd) !rxxfreed correspond to : !som[(exp(-BEi/Z)).{(S-,ij)(S+,ji)+(S-,ij)(S+,ji)}/2{0.5Jf(wI-wij,td,tn)+3Jf(wI+wij,td,tn)}] spijsmji= mdens(ni)*(szpld(ni,nj)*sznegd(nj,ni)+ & sznegd(ni,nj)*szpld(nj,ni))/2 dxxfreed=spijsmji*(0.5*rfie+3.*rfien) rxxfreed=rxxfreed + dxxfreed !rzxfreed correspond som[(exp(-BEi/Z)).(Sz,ij).(Sz,ji).{Jf(wI-wij,td,tn)+Jf(wI+wij,td,tn)}/2] szijszji=mdens(ni)*szzd(ni,nj)*szzd(nj,ni) dzxfreed=szijszji*(rfie+rfien)/2. rzxfreed=rzxfreed + dzxfreed !rzzfreed correspond to: som[(exp(-BEi/Z)).(Sz,ij).(Sz,ji).{2*Jf(wI,td,tn)}] dzzfreed=szijszji*2.*rfe rzzfreed=rzzfreed + dzzfreed !rxzfreed correspond to :som[(exp(-BEi/Z)).{(S-,ij)(S+,ji)+(S-,ij)(S+,ji)}/2{1.5Jf(wI,td,tn)}] dxzfreed=spijsmji*(rfe) rxzfreed=rxzfreed+dxzfreed end do rfei=rfreed(cmt,omegaii,taun1,tcd) rfez=rfreed(cmt,0.,taun1,tcd) !rayan1 correspond to the term -3sz*sz*(Jf(wI,td,tn) -Ja(wI,tD)) !rayan2 correspond to the term -2sz*sz*(Jf(0,td,tn)+3/5Jf(wI,td,tn)-Ja(0,tD)-3/4Ja(wI,tD)) rayan1=rayan1-3.0*szm*szm*(rfei-raya)*mdens(ni) rayan2=rayan2-2.0*szm*szm* & (rfez+0.75*rfei-rayz-0.75*raya)*mdens(ni) end do rayan1=rayan1*mdensangl(nagl) rayan2=rayan2*mdensangl(nagl) r2(2)= r2(2) + rayan2 r1(2)= r1(2) + rayan1 r1freed = mdensangl(nagl)*(rxxfreed + 3*rzxfreed) r2freed = mdensangl(nagl)*(0.5*rxxfreed + 1.5*rzxfreed + & rzzfreed + 1.5*rxzfreed) r1(1)=r1(1)+r1freed r2(1)=r2(1)+r2freed if (sometang.ne.0) then end if 6000 continue r1t=r1(1)+r1(2) r2t=r2(1)+r2(2) WRITE (*,8001) 1000*CHBO,r1t/sometang,r2t/sometang WRITE (3,8001) 1000*CHBO,r1t/sometang,r2t/sometang out=r1t/sometang write (4,8002) 1000*CHBO,r1t/sometang !5000 continue 7000 continue CLOSE(3) CLOSE(4) ! write (*,*)tcd,taun 8001 format(f9.4,f9.4,f9.4) 8002 format(f9.4,f9.4) return END FUNCTION RFREED(CMT,OMEGA,TAUN,TCD) TETA=ATAN(OMEGA*TAUN) RHO=TCD*SQRT(OMEGA*OMEGA+1./(TAUN*TAUN)) CO1=COS(TETA/2.) CO2=COS(TETA) CO3=COS(TETA*3./2.) SI1=SIN(TETA/2.) SI2=SIN(TETA) SI3=SIN(TETA*3./2.) RHO2=SQRT(RHO) RHO32=RHO2*RHO2*RHO2 xNUM=1.+RHO/4.+ RHO2*5.*CO1/4.+(4.*RHO*CO2)/9.+RHO32* & (CO3+CO1)/9.+(CO2*RHO*RHO)/36. DEN1=1.+RHO2*CO1+4./9.*RHO*CO2+RHO32*CO3/9. DEN2=RHO2*SI1+4./9.*RHO*SI2+SI3/9.*RHO32 RFREED=3.*cmt*xNUM/(DEN1*DEN1+DEN2*DEN2) END FUNCTION RFREED FUNCTION RAYANT(CMT,OMEGA,TCD) Z=SQRT(2.*TCD*OMEGA) XJSUP=1.+5.*Z/8.+Z*Z/8. XJINF=1.+Z+Z*Z/2.+Z*Z*Z/6.+4.*Z*Z*Z*Z/81.+Z*Z*Z*Z*Z/81.+ & Z*Z*Z*Z*Z*Z/648. RAYANT=3.*cmt*XJSUP/XJINF ! write(3,*) 'ayant w r',omega,rayant END FUNCTION RAYANT SUBROUTINE FUNCsuper(P,FUNC,NMX,NV) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(NV) DIMENSION XI2(21,21) real*4 chbo,out COMMON /SET/SET COMMON /PPAR/ P2(21) COMMON /RK10/ SPIN, SI COMMON /WATER/ ACQ COMMON /STAMPA/INDEXSTAMPA COMMON /STEPGAMMA/ STEPGAMMA COMMON /B4/ NVMEM,NPT(10),NPTOT COMMON /B31/ TPUNO(500) COMMON /B32/ PP(500),Z(500) COMMON /C1M/DM(10),DDM(10),CONCM(10) COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10), & TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),GXM(10), & GYM(10),GZM(10) COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10), & TAUMM(10,10) COMMON /MOLFRAZM/ AMOLFRAM(10) COMMON /CONTATM/ ACONTM(10) COMMON /INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10) & ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10) COMMON /C1/D,DD,CONC COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4 COMMON /GTENSOR/ GX,GY,GZ COMMON /TAU1/ TAUS0 COMMON /TAU/ TAUR,TAUV,TAUM COMMON /MOLFRAZ/ AMOLFRA COMMON /CONTAT/ ACONT COMMON /GAMMAH/ GAMMAI COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS COMMON /CICLE/ NVEST COMMON /TMAT/ TMAT(500,10),TMATCONT(500,10),TMATCROSS(500,10) COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10), & B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10), & B16(10),B17(10),B18(10),B19(10),B20(10),B21(10) COMMON/TEMPERATURE/ TEMP(10) DIMENSION TS1(10),TS2(10),TR1(10),TR2(10),TV1(10),TV2(10) COMMON /TMSTART/ TM11(10),TM21(10) common /super/ b,DIF,xMSAT,TAUN,HANI,TEMP1 common /super2/ omega0,ppar,qpar real*4 omega0,ppar,qpar real*4 b,DIF,xMSAT,TAUN,HANI,TEMP1 C SET PARAMETERS FB=0. ind=1 IF(B1(1).EQ.1)then TAUN=P(ind) ind=ind+1 endif if(b13(1).eq.1)then b=p(ind) ind=ind+1 endif if(b13(1).eq.2)then ppar=p(ind) ind=ind+1 endif if(b13(1).eq.3)then qpar=p(ind) ind=ind+1 endif if(b13(1).eq.4)then b=p(ind) ind=ind+1 ppar=p(ind) ind=ind+1 qpar=p(ind) ind=ind+1 endif if(b13(1).eq.5)then ppar=p(ind) ind=ind+1 qpar=p(ind) ind=ind+1 endif if(b14(1).eq.1)then dif=p(ind) ind=ind+1 endif if(b14(1).eq.2)then xmsat=p(ind) ind=ind+1 endif if(b14(1).eq.3)then hani=p(ind) ind=ind+1 endif if(b14(1).eq.4)then dif=p(ind) ind=ind+1 xmsat=p(ind) ind=ind+1 hani=p(ind) ind=ind+1 endif if(b14(1).eq.5)then xmsat=p(ind) ind=ind+1 hani=p(ind) ind=ind+1 endif if(b14(1).eq.6)then dif=p(ind) ind=ind+1 hani=p(ind) ind=ind+1 endif if(b14(1).eq.7)then dif=p(ind) ind=ind+1 xmsat=p(ind) ind=ind+1 endif omega0=hani C PARAMETERS************************************************************** do i=1,npt(1) chbo=z(i)/1000. ! write(6,*)ppar,qpar CALL provag2(chbo,out) FB=((PP(I)-out)**2)/abs(out)+FB IF(INDEXSTAMPA.EQ.1)WRITE(44,'(2X,2(E10.4))')Z(I),out end do ! IF(INDEXSTAMPA.EQ.1)then ! do i=1,nv ! write(6,*)p(i) ! end do ! endif ! IF(INDEXSTAMPA.EQ.1)WRITE(6,*)' ERROR', ! & '=',FB FUNC=fb RETURN END SUBROUTINE FUNCsuper2(P,FUNC,NMX,NV) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(NV) DIMENSION XI2(21,21) real*4 chbo,out COMMON /SET/SET COMMON /PPAR/ P2(21) COMMON /RK10/ SPIN, SI COMMON /WATER/ ACQ COMMON /STAMPA/INDEXSTAMPA COMMON /STEPGAMMA/ STEPGAMMA COMMON /B4/ NVMEM,NPT(10),NPTOT COMMON /B31/ TPUNO(500) COMMON /B32/ PP(500),Z(500) COMMON /C1M/DM(10),DDM(10),CONCM(10) COMMON /IPERFM/AZM(10),AYM(10),AXM(10),THETAM(10),RKM(10), & TAUCM(10),DPARAM(10),EPARAM(10),PHIM(10),S4M(10),GXM(10), & GYM(10),GZM(10) COMMON /TAU1M/ TAUS0M(10,10),TAURM(10,10),TAUVM(10,10), & TAUMM(10,10) COMMON /MOLFRAZM/ AMOLFRAM(10) COMMON /CONTATM/ ACONTM(10) COMMON /INDA/ INDDPARA(10),INDAPAR(10),INDAPER(10),INDEPARA(10) & ,INDED(10),INDAMOLFRA(10),INDS4(10),INDAPER2(10) COMMON /C1/D,DD,CONC COMMON /IPERF/AZ,AY,AX,THETA,RK,TAUC,DPARA,EPARA,PHI,S4 COMMON /GTENSOR/ GX,GY,GZ COMMON /TAU1/ TAUS0 COMMON /TAU/ TAUR,TAUV,TAUM COMMON /MOLFRAZ/ AMOLFRA COMMON /CONTAT/ ACONT COMMON /GAMMAH/ GAMMAI COMMON /TM/ TMUNO,TMUNOCONT,TMUNOCROSS COMMON /CICLE/ NVEST COMMON /TMAT/ TMAT(500,10),TMATCONT(500,10),TMATCROSS(500,10) COMMON /BPARA/B1(10),B2(10),B3(10),B4(10),B5(10),B6(10),B7(10), & B8(10),B9(10),B10(10),B11(10),B12(10),B13(10),B14(10),B15(10), & B16(10),B17(10),B18(10),B19(10),B20(10),B21(10) COMMON/TEMPERATURE/ TEMP(10) DIMENSION TS1(10),TS2(10),TR1(10),TR2(10),TV1(10),TV2(10) COMMON /TMSTART/ TM11(10),TM21(10) common /super/ b,DIF,xMSAT,TAUN,HANI,TEMP1 common /super2/ omega0,ppar,qpar real*4 omega0,ppar,qpar real*4 b,DIF,xMSAT,TAUN,HANI,TEMP1 C SET PARAMETERS FB=0. ind=1 IF(B1(1).EQ.1)then TAUN=P(ind) ind=ind+1 endif if(b13(1).eq.1)then b=p(ind) ind=ind+1 endif if(b14(1).eq.1)then dif=p(ind) ind=ind+1 endif if(b14(1).eq.2)then xmsat=p(ind) ind=ind+1 endif if(b14(1).eq.3)then hani=p(ind) ind=ind+1 endif if(b14(1).eq.4)then dif=p(ind) ind=ind+1 xmsat=p(ind) ind=ind+1 hani=p(ind) ind=ind+1 endif if(b14(1).eq.5)then xmsat=p(ind) ind=ind+1 hani=p(ind) ind=ind+1 endif if(b14(1).eq.6)then dif=p(ind) ind=ind+1 hani=p(ind) ind=ind+1 endif if(b14(1).eq.7)then dif=p(ind) ind=ind+1 xmsat=p(ind) ind=ind+1 endif C PARAMETERS************************************************************** do i=1,npt(1) chbo=z(i)/1000. ! write(6,*) b,DIF,xMSAT,TAUN,HANI,TEMP1,chbo CALL provag (chbo,out) FB=((PP(I)-out)**2)/abs(out)+FB IF(INDEXSTAMPA.EQ.1)WRITE(44,'(2X,2(E10.4))')Z(I),out end do IF(INDEXSTAMPA.EQ.0)WRITE(6,*)' ERROR', & '=',FB FUNC=fb RETURN end subroutine provag2(chbo,out) real*4 out,chbo real*4 b,DIF,xMSAT,TAUN,HANI,TEMP real*4 omega0,ppar,qpar real*4 l,x,omega,omegas,omegai,jf1,jf2,ja,jf3 real*4 tcd,xmaspar,sisatu,nelec,spinmul,xmsatu,concpar,cmt real*4 chbot common /super/ b,DIF,xMSAT,TAUN,HANI,TEMP common /super2/ omega0,ppar,qpar ! write(6,*) b,DIF,xMSAT,TAUN,HANI,TEMP1,chbo,ppar,qpar,omega0 hskghz=0.06626/1.3181 ! spinmul ratio between real and fictif spin TCD=(B*B)/dif xMASPAR=3.1416*4*B*B*B*5.18/3 sisatu=xmsat*xmaspar/1000 nelec=sisatu/9.274e-24 spinmul=nelec/2 xMSATU=xMSAT*xMASPAR CONCPAR=232E-3/(3*xMASPAR) CMT=1.7765e+5*CONCPAR*xMSATU*xMSATU/(b*dif) chbot=chbo*4.7/0.2 !campo in T !23.81 x=sisatu*chbot/(1.3806e-23*temp) !k cothx=1./tanh(x) l=cothx-1./x ! write(6,*)x,l omegas=chbo*1000000000.*6.28*658 omegai=chbo*1000000000.*6.28 if(omegas.ge.omega0)then omega=(omegas**(1./4)-omega0**(1./4))**4 else omega=0 endif TCD=(B*B)/dif ! write(6,*)cmt,omega,taun,tcd jf1=rfreed(cmt,omega,taun,tcd) jf2=rfreed(cmt,omegai,taun,tcd) jf3=rfreed(cmt,0,taun,tcd) ja=rayant(cmt,omegai,tcd) jf1=jf1/3 jf2=jf2/3 jf3=jf3/3 ja=ja/3 ! write(6,*)ppar,qpar,jf3/3,cmt*(7./3*(ppar+qpar)+1)*jf3/3 out=cmt*(7*ppar*l/x*jf1+(7*qpar*l/x+3*(ppar+qpar)* & (1.-l*l-2*l/x))*jf2+3*l*l*ja) out2=cmt/2*(13*ppar*l/x*jf1+7*qpar*l/x*jf2+6*qpar*l/x*jf3+ & (ppar+qpar)*(1.-l*l-2*l/x)*(3*jf2+4*jf3)+l*l*(3*ja+4.)) ! write(6,*)out return end