SUBROUTINE H2H2 C C FARRAR-LEE H2-H2 MSV POTENTIAL PLUS C ANISOTROPIC TERMS USED BY ZARUR AND RABITZ. C C ENTRY VINIT(I,RM,EPSIL) INITIALIZES ROUTINE FOR LAM(I) C SETS RM, EPSIL C ENTRY VSTAR (I,R,V) GIVES POTENTIAL FOR LAM(I) AT R. C ENTRY VSTAR1(I,R,V) GIVES DERIVATIVE FOR LAM(I) AT R. C ENTRY VSTAR2(I,R,V) GIVES 2ND DERIV. FOR LAM(I) AT R. IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION B(4),A(4,4) DIMENSION CONST(4) C C DEFINE STATEMENT FUNCTIONS E(R)=EXP(BETA*(1D0-R)) P4(C1,C2,C3,C4,R)=C1+R*(C2+R*(C3+R*C4)) P3(C1,C2,C3,R)=C1+R*(C2+R*C3) P2(C1,C2,R)=C1+R*C2 C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ENTRY VINIT(I,RM,EPSIL) C C SET CONSTANT FACTORS . . . C RM=3.49D0 EPSIL=24.17D0 RMSAVE=RM CONST(1)=(4D0*3.1415926D0)**( 1.5D0) CONST(2)=(.14D0/5D0)*CONST(1) CONST(3)=CONST(2) CONST(4)=54836D0 / EPSIL C IF (I.EQ.4) GO TO 1400 IF (I.LT.1 .OR. I.GT.4) GO TO 9999 C WRITE(6,601) 601 FORMAT('0',10X,'FARRAR-LEE MSV H2-H2 POTENTIAL / RM=3.49,', 1 ' EPSIL=24.17 SET INTERNALLY.') IF (I.GT.1) WRITE(6,606) 606 FORMAT(16X,'ANISOTROPIC TERM IS .14/5.0 TIMES ISOTROPIC TERM') C C THE FOLLOWING PARAMETERS MUST BE SPECIFIED INTERNALLY. C R1,R2 IN RM . C6,C8 IN (1/CM )/ANG**N BETA C C SCALE PARAMETERS WITH RM, EPS. C R1=1.0754D0 R2=1.4D0 C6=( 57913.D0 /EPSIL)/ RM**6 C8=( 156263.D0/EPSIL)/ RM**8 BETA=6.5D0 C WRITE(6,602) R1,BETA,R2,C6,C8 602 FORMAT(20X,'FOR R LESS THAN',F8.4,', MORSE BETA =',F9.4/ 1 20X,'FOR R GREATER THAN',F8.4,', VAN DER WAAL C6, C8 =', 2 2F10.5,' (EPSIL/RM**N).') C C SET-UP SPLINE COEFFS. B(J,I), J=1,4 FOR LAM(I). C ** MAR 1979 CHANGED TO INCORPORATE PRECOMPUTED SPLINE COEFFS. C B(1)=-.6829 4444 3121 D1 B(2)= .7211 3009 9412 D1 B(3)=-.7689 9426 2915 D0 B(4)=-.7125 2209 5040 D0 C 1100 WRITE(6,605) (B(J ),J=1,4) 605 FORMAT(20X,'CUBIC SPLINE COEFFICIENTS, TO INCREASING POWERS OF R, 1ARE AS FOLLOWS'/25X,4D16.8) RETURN C 1400 WRITE(6,607) 607 FORMAT('0',10X,'QUADRUPOLE-QUADRUPOLE LONG-RANGE W/ Q=.662 BUC 1KINGHAM.') AC=.0687D0 RETURN C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ENTRY VSTAR (I,R,V) IF (I.EQ.4) GO TO 2400 IF (R.LE.R1) GO TO 2100 IF (R.GE.R2) GO TO 2200 V=P4(B(1),B(2),B(3),B(4),R) GO TO 5000 2100 T1=E(R) V=T1*(T1-2D0) GO TO 5000 2200 RSQ=1D0/(R*R) R6=RSQ*RSQ*RSQ V=- R6*(C6+C8*RSQ) GO TO 5000 99. 2400 RR=R*RMSAVE R12=R**(-12) RRT=RR+AC*R12 V=RRT**(-5) GO TO 5000 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ENTRY VSTAR1(I,R,V) IF (I.EQ.4) GO TO 3400 IF (R.LE.R1) GO TO 3100 IF (R.GE.R2) GO TO 3200 V=P3(B(2),2D0*B(3),3D0*B(4),R) GO TO 5000 3100 T1=E(R) V=-2D0*BETA* T1*(T1-1D0) GO TO 5000 3200 RSQ=1D0/(R*R) R6=RSQ*RSQ*RSQ/R V= R6*(6D0*C6+8D0*C8*RSQ) GO TO 5000 3400 RR=R*RMSAVE R12=R**(-12) RRT=RR+AC*R12 V=-5D0*RRT**(-6)*(1D0-12D0*AC*R12/RR) GO TO 5000 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C ENTRY VSTAR2(I,R,V) IF (I.EQ.4) GO TO 4400 IF (R.LE.R1) GO TO 4100 IF (R.GE.R2) GO TO 4200 V=P2(2D0*B(3),6D0*B(4),R) GO TO 5000 4100 T1=E(R) V=T1*BETA*BETA* (4D0*T1-2D0) GO TO 5000 4200 RSQ=1D0/(R*R) R6=RSQ*RSQ R6=R6*R6 V=- R6*(42D0*C6+56D0*C8*RSQ) GO TO 5000 4400 RR=R*RMSAVE R12=R**(-12) RRT=RR+AC*R12 XXX=1D0-12D0*AC*R12/RR V=-5D0*RRT**(-6)*(156D0*AC*R12/(RR*RR)-6D0*XXX*XXX/RRT) GO TO 5000 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C UNIFIED RETURN POINT 5000 V=V*CONST(I) RETURN C C 9999 WRITE(6,699) I 699 FORMAT('0 * * * ERROR. POTENTIAL NOT DEFINED FOR SYMMETRY=',I10) STOP END SUBROUTINE VRTP(IDERIV,R,VV) C C FARRAR-LEE H2-H2 MSV POTENTIAL PLUS C ANISOTROPIC TERMS USED BY ZARUR AND RABITZ. C COMMON /ANGLES/ FOR VERSION 14 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE DIMENSION B(4),A(4,4) DIMENSION CONST(4),VI(4),FACT(4) C COMMON/ANGLES/COSANG(7),FACTOR,IHOMO,ICNSYM,IHOMO2,ICNSY2 C C DEFINE STATEMENT FUNCTIONS E(R)=EXP(BETA*(1D0-R)) P4(C1,C2,C3,C4,R)=C1+R*(C2+R*(C3+R*C4)) P3(C1,C2,C3,R)=C1+R*(C2+R*C3) P2(C1,C2,R)=C1+R*C2 C C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IF (IDERIV.EQ.-1) THEN RM=3.49D0 RMSAVE=RM EPSIL=24.17D0 C SET UP CONSTANTS CONST(1)=(4D0*3.1415926D0)**( 1.5D0) CONST(2)=(.14D0/5D0)*CONST(1) CONST(3)=CONST(2) CONST(4)=54836D0 / EPSIL R1=1.0754D0 R2=1.4D0 C6=( 57913.D0 /EPSIL)/ RM**6 C8=( 156263.D0/EPSIL)/ RM**8 BETA=6.5D0 AC=.0687D0 C SET-UP SPLINE COEFFS. B(J,I), J=1,4 FOR LAM(I). B(1)=-.6829 4444 3121 D1 B(2)= .7211 3009 9412 D1 B(3)=-.7689 9426 2915 D0 B(4)=-.7125 2209 5040 D0 WRITE(6,601) RM,EPSIL 601 FORMAT(10X,'FARRAR-LEE MSV H2-H2 POTENTIAL'/10X,'RM=',F8.3/ 1 10X,'EPSIL=',F12.2) WRITE(6,606) 606 FORMAT(10X,'ANISOTROPIC TERM IS .14/5.0 TIMES ISOTROPIC TERM') WRITE(6,602) R1,BETA,R2,C6,C8 602 FORMAT(10X,'FOR R LESS THAN',F8.4,', MORSE BETA =',F9.4/ 1 10X,'FOR R GREATER THAN',F8.4,', VAN DER WAAL C6, C8 =', 2 2F10.5,' (EPSIL/RM**N).') WRITE(6,607) 607 FORMAT(10X,'QUADRUPOLE-QUADRUPOLE LONG-RANGE W/ Q=.662 BU', 1 'CKINGHAM.') WRITE(6,605) (B(J),J=1,4) 605 FORMAT(10X,'CUBIC SPLINE COEFFICIENTS, TO INCREASING ', 1 'POWERS OF R, ARE AS FOLLOWS'/15X,4D16.8) PI=ACOS(-1.D0) PIF=4.D0*PI FACT(1)=1.D0/SQRT(PIF*PIF*PIF) FACT(2)=FACT(1)*2.5D0 FACT(3)=FACT(1) FACT(4)=FACT(1)*(45.D0/4.D0)*SQRT(1.D0/7.D1) C IHOMO=2 ICNSYM=2 C RETURN RM EPSIL R=RM VV=EPSIL RETURN C ELSEIF (IDERIV.EQ.0) THEN C IF (R.LE.R1) GO TO 2100 IF (R.GE.R2) GO TO 2200 V=P4(B(1),B(2),B(3),B(4),R) GO TO 5000 2100 T1=E(R) V=T1*(T1-2D0) GO TO 5000 2200 RSQ=1D0/(R*R) R6=RSQ*RSQ*RSQ V=-R6*(C6+C8*RSQ) 5000 DO 5001 I=1,3 5001 VI(I)=V*CONST(I) C QUADRUPOLE - QUADRUPOLE TERM RR=R*RMSAVE R12=R**(-12) RRT=RR+AC*R12 V=RRT**(-5) VI(4)=V*CONST(4) C CT1=COSANG(1) C ST1=SQRT(1.D0-CT1*CT1) C CT2=COSANG(2) C ST2=SQRT(1.D0-CT2*CT2) C P21=3.D0*CT1*CT1-1.D0 C P22=3.D0*CT2*CT2-1.D0 C ADD THE PIECES WITH COSANG() VALUES - YRR GIVES ROTOR FUNCTIONS VV=VI(1)*YRR(0,0,0,COSANG(1),COSANG(2),COSANG(3)) 1 +VI(2)*YRR(2,0,2,COSANG(1),COSANG(2),COSANG(3)) 2 +VI(3)*YRR(0,2,2,COSANG(1),COSANG(2),COSANG(3)) 1 +VI(4)*YRR(2,2,4,COSANG(1),COSANG(2),COSANG(3)) RETURN ELSE C WRITE(6,*) ' *** VRTP (H2-H2). NO SUPPORT FOR IDERIV',IDERIV STOP ENDIF END