subroutine HHO_ESA_init() c V. Engel and R. Schinke, J.C.P vol. 88, 129 (1988) return end c================================== subroutine HHO_ESA_pes(xx,pe) implicit none real*8 xx(3),x(3),pot,pe c transform coordinates from HHO to HOH x(1)=xx(2) x(2)=xx(1) x(3)=xx(3) x(1)=x(1)*0.5292d0 x(2)=x(2)*0.5292d0 x(3)=x(3)*0.5292d0 pe=pot(x)/27.2116d0 x(1)=x(1)/0.5292d0 x(2)=x(2)/0.5292d0 x(3)=x(3)/0.5292d0 return end c--------- FUNCTION POTI(r) C C H2O - FITPOTENTIAL C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3),RE(3),ALPHA(3),C(50),PHI(100) C DATA DOH,BETOH,ROH/4.621D0 , 2.294D0 , 0.971D0/ data norder/50/ data re/1.0000d0,1.0000d0,1.0000d0/ data alpha/1.2500d0,0.7500d0,1.2500d0/ data gamma/0.25d0/ data c/ -.93562626d0, -.33992392d1, -.37923375d1, $ .35229926d1, .20508566d1, .33290770d0, $ .64682997d0, -.51447712d2, -.49046646d1, $ .24435891d2, .76337272d0, .10582765d2, $ .73669551d1, .64329848d2, .89923222d1, $ -.89690335d2, .64529244d2, -.53894246d2, $ -.59093483d0, -.12138697d2, .65697698d1, $ .59787268d2, -.29842680d2, .30342123d2, $ -.29963985d2, .24344415d2, -.62823352d2, $ .18190455d3, .10358976d2, -.39393237d2, $ .39698652d1, .12165069d2, .51422343d1, $ -.47491777d1, -.15323205d1, .92191236d1, $ -.74335345d1, .15738989d3, .33299372d1, $ -.35652609d1, -.15416913d3, -.17472317d1, $ .52101774d2, .16650884d3, -.23987576d2, $ -.69990914d2, .19884456d2, .40429481d2, $ -.97198548d1, .18539893d1/ c--------------------------- C RMORSE(X)=DOH*(DEXP(-BETOH*(X-ROH))-1.)**2 - DOH XSI(X,Y)=(1.-DEXP(-gamma*X**2)) * (1.-DEXP(-gamma*Y**2)) C C EINLESEN DER DATEN C c READ(22,100)NORDER,(RE(I),I=1,3),(ALPHA(I),I=1,3),gamma c *,(C(I),I=1,NORDER) c 100 format(i3,/,3f12.4,/,3f12.4,/,f12.4,50(/,e16.8)) write(6,*) 'Nach Einlesen der Daten norder:',norder DUMMY=RINTI(R) RETURN C ENTRY POT(r) C RMAX=DMAX1(R(1),R(3))/0.5292 IF(RMAX.GE.5.5)GOTO 10 CALL FUNC(R,PHI,NORDER,RE,ALPHA) POT=0.0 DO 1 I=1,NORDER POT=POT+C(I)*PHI(I) 1 continue POT=POT+RMORSE(R(1))*XSI(R(2),R(3)) * +RMORSE(R(3))*XSI(R(2),R(1)) IF(RMAX.LE.3.5)RETURN VINT=RINT(R) COSCOS=DCOS(1.5707963*(RMAX-3.5)) GW=0.5+0.5*COSCOS FW=0.5-0.5*COSCOS POT=POT*GW+VINT*FW RETURN C 10 CONTINUE POT=RINT(R) RETURN END c FUNCTION RINTI(R) C C H2O INTERPOLATIONS POTENTIAL C IMPLICIT REAL*8 (A-H,O-Z) REAL RLGDRE C DIMENSION R(3),GAMMA(5),A(5,5),B(5,6),P(5),X(6) C DATA D0,ALPHA0,R0/4.621,2.294,0.971/ DATA GAMMA/45.,75.,104.52,135.,180./ DATA B/-0.69629,-0.48884,-0.45792,-0.44409,-0.41761, * 1.40886,1.39153,1.39691,1.46819,1.54303, * 0.12042,0.14350,0.17865,0.19215,0.19151, * 0.19193,0.15280,0.13400,0.09861,0.05140, * 0.62717,1.38163,1.43055,1.36016,1.22214, * 0.31021,0.12674,0.18001,0.38321,0.51142/ C C C BESTIMMUNG DER LEGENDRE ENTWICKLUNGS KOEFF. C DO 1 I=1,5 DO 1 J=1,5 1 A(I,J)=RLGDRE(J-1,COS(1.7453292E-2*SNGL(GAMMA(I)))) CALL MAT56(A,5,B,6,DET) RETURN C ENTRY RINT(R) C RMAX=DMAX1(R(1),R(3)) RMAXA0=RMAX/0.5292 - 3.5 RMIN=DMIN1(R(1),R(3)) COSGAM=-(R(2)**2-R(1)**2-R(3)**2)/(2*R(1)*R(3)) DO 2 I=1,5 2 P(I)=RLGDRE(I-1,SNGL(COSGAM)) DO 3 J=1,6 X(J)=0.0 DO 3 I=1,5 X(J)=X(J)+B(I,J)*P(I) 3 CONTINUE D=D0+X(1)*DEXP(-X(2)*RMAXA0-X(3)*RMAXA0**2) ALPHA=ALPHA0+X(4)*DEXP(-X(5)*RMAXA0-X(6)*RMAXA0**2) RINT=D*(DEXP(-ALPHA*(RMIN-R0))-1)**2-D RETURN END c FUNCTION RLGDRE(N,X) C CALCULATES LEGENDRE POLYNOMIALS USING RECURRENCE ABRAM./STEGUN C PAGE 782 RLGDRE=1 R2=RLGDRE IF(N.EQ.0)RETURN RLGDRE=X R1=RLGDRE IF(N.EQ.1)RETURN DO 1 I=2,N RXR1=X*R1 RLGDRE=RXR1+RXR1-R2+(R2-RXR1)/I R2=R1 1 R1=RLGDRE RETURN END c c c SUBROUTINE MAT56(A,N,B,M,DETERM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION INDEX(5,2),IPIVOT(5),A(5,5),B(5,6),PIVOT(5) EQUIVALENCE (IROW,JROW), (ICOLUM,JCOLUM), (AMAX, T, SWAP) C INITIALIZATION 10 DETERM=1.0D0 15 DO 20 J=1,N 20 IPIVOT(J)=0 30 DO 550 I=1,N C SEARCH FOR PIVOT ELEMENT 40 AMAX=0.0D0 45 DO 105 J=1,N 50 IF (IPIVOT(J)-1) 60, 105, 60 60 DO 100 K=1,N 70 IF (IPIVOT(K)-1) 80, 100, 740 80 IF( DABS(AMAX)- DABS(A(J,K)))85,100,100 85 IROW=J 90 ICOLUM=K 95 AMAX=A(J,K) 100 CONTINUE 105 CONTINUE 110 IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1 C INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 130 IF (IROW-ICOLUM) 140, 260, 140 140 DETERM=-DETERM 150 DO 200 L=1,N 160 SWAP=A(IROW,L) 170 A(IROW,L)=A(ICOLUM,L) 200 A(ICOLUM,L)=SWAP 205 IF(M) 260, 260, 210 210 DO 250 L=1, M 220 SWAP=B(IROW,L) 230 B(IROW,L)=B(ICOLUM,L) 250 B(ICOLUM,L)=SWAP 260 INDEX(I,1)=IROW 270 INDEX(I,2)=ICOLUM 310 PIVOT(I)=A(ICOLUM,ICOLUM) 320 DETERM=DETERM*PIVOT(I) C DIVIDE PIVOT ROW BY PIVOT ELEMENT 330 A(ICOLUM,ICOLUM)=1.0D0 340 DO 350 L=1,N 350 A(ICOLUM,L)=A(ICOLUM,L)/PIVOT(I) 355 IF(M) 380, 380, 360 360 DO 370 L=1,M 370 B(ICOLUM,L)=B(ICOLUM,L)/PIVOT(I) C REDUCE NON-PIVOT ROWS 380 DO 550 L1=1,N 390 IF(L1-ICOLUM) 400, 550, 400 400 T=A(L1,ICOLUM) 420 A(L1,ICOLUM)=0.0D0 430 DO 450 L=1,N 450 A(L1,L)=A(L1,L)-A(ICOLUM,L)*T 455 IF(M) 550, 550, 460 460 DO 500 L=1,M 500 B(L1,L)=B(L1,L)-B(ICOLUM,L)*T 550 CONTINUE C INTERCHANGE COLUMNS 600 DO 710 I=1,N 610 L=N+1-I 620 IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630 630 JROW=INDEX(L,1) 640 JCOLUM=INDEX(L,2) 650 DO 705 K=1,N 660 SWAP=A(K,JROW) 670 A(K,JROW)=A(K,JCOLUM) 700 A(K,JCOLUM)=SWAP 705 CONTINUE 710 CONTINUE 740 RETURN END c c c SUBROUTINE FUNC(R,PHI,NORDER,X,ALPHA) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3),S(3),S2(3),S3(3),S4(3),S5(3) DIMENSION PHI(100),S6(3) DIMENSION ALPHA(3),X(3) COMMON BETA(3) C FAK=1.0 DO 1 I=1,3 S(I)=R(I)-X(I) S2(I)=S(I)**2 S3(I)=S2(I)*S(I) S4(I)=S2(I)*S2(I) S5(I)=S3(I)*S2(I) S6(I)=S3(I)**2 FAK=FAK*(1.-DTANH(ALPHA(I)*S(I))) 1 CONTINUE C C PHI(1)=1. C PHI(2)=S(1)+S(3) PHI(3)=S(2) C PHI(4)=S2(1)+S2(3) PHI(5)=S2(2) PHI(6)=S(2)*(S(1)+S(3)) PHI(7)=S(1)*S(3) C PHI(8)=S3(1)+S3(3) PHI(9)=S3(2) PHI(10)=S(1)*S2(3)+S2(1)*S(3) PHI(11)=S2(2)*(S(1)+S(3)) PHI(12)=S(2)*(S2(1)+S2(3)) PHI(13)=S(1)*S(2)*S(3) C PHI(14)=S4(1)+S4(3) PHI(15)=S4(2) PHI(16)=S3(1)*S(3)+S(1)*S3(3) PHI(17)=S2(1)*S2(3) PHI(18)=S(2)*(S3(1)+S3(3)) PHI(19)=S2(2)*(S2(1)+S2(3)) PHI(20)=S3(2)*(S(1)+S(3)) PHI(21)=S2(2)*S(1)*S(3) PHI(22)=S(2)*(S2(1)*S(3)+S(1)*S2(3)) C C PHI(23)=S5(1)+S5(3) PHI(24)=S4(1)*S(3)+S(1)*S4(3) PHI(25)=S3(1)*S2(3)+S2(1)*S3(3) PHI(26)=S(2)*(S4(1)+S4(3)) PHI(27)=S(2)*(S3(1)*S(3)+S(1)*S3(3)) PHI(28)=S(2)*S2(1)*S2(3) PHI(29)=S2(2)*(S3(1)+S3(3)) PHI(30)=S2(2)*(S2(1)*S(3)+S(1)*S2(3)) PHI(31)=S3(2)*(S2(1)+S2(3)) PHI(32)=S3(2)*S(1)*S(3) PHI(33)=S4(2)*(S(1)+S(3)) PHI(34)=S5(2) C PHI(35)=S6(1)+S6(3) PHI(36)=S5(1)*S(3)+S(1)*S5(3) PHI(37)=S4(1)*S2(3)+S2(1)*S4(3) PHI(38)=S3(1)*S3(3) PHI(39)=S(2)*(S5(1)+S5(3)) PHI(40)=S(2)*(S4(1)*S(3)+S(1)*S4(3)) PHI(41)=S(2)*(S3(1)*S2(3)+S2(1)*S3(3)) PHI(42)=S2(2)*(S4(1)+S4(3)) PHI(43)=S2(2)*(S3(1)*S(3)+S(1)*S3(3)) PHI(44)=S2(1)*S2(2)*S2(3) PHI(45)=S3(2)*(S3(1)+S3(3)) PHI(46)=S3(2)*(S2(1)*S(3)+S(1)*S2(3)) PHI(47)=S4(2)*(S2(1)+S2(3)) PHI(48)=S4(2)*S(1)*S(3) PHI(49)=S5(2)*(S(1)+S(3)) PHI(50)=S6(2) C C DO 2 I=1,NORDER PHI(I)=FAK*PHI(I) 2 CONTINUE C RETURN END C C FUNCTION SPHARM(L,M,X,Y) C CALCULATES SPAERICAL HARMONICS EDMONDS 2.5.29 COMPLEX SPHARM REAL*8 FAK,ASLEG IF(M.NE.0)GOTO 1 SPHARM=0.2820947917*SQRT(L+L+1.)*RLGDRE(L,COS(X)) **CMPLX(1.,0.) RETURN 1 CONTINUE MABS=IABS(M) IF(MABS.LE.L)GOTO 3 spharm=0.0 return 3 CONTINUE IO=MABS+MABS FAK=1.D0 DO 2 I=1,IO 2 FAK=FAK*SQRT(FLOAT(L-MABS+I)) YM=MABS*Y SPHAR=0.2820947917*SQRT(L+L+1.)/SNGL(FAK)*SNGL(ASLEG(L,MABS,X))* *(-1.)**MABS SPHARM=SPHAR*CMPLX(COS(YM),SIN(YM)) IF(M.GE.0)RETURN SPHARM=(-1)**M*CONJG(SPHARM) RETURN END C C C FUNCTION ASLEG(L,M,BETA) REAL*8 ASLEG,FAK,PLM1,PLM2 X=COS(BETA) MABS=IABS(M) IF(MABS.GT.0)GOTO 1 ASLEG=RLGDRE(L,X) RETURN 1 FAK=1.D0 DO 2 I=1,MABS 2 FAK=FAK*(MABS+I) ARG=0.5*BETA PLM2=FAK*(COS(ARG)*SIN(ARG))**MABS IF(L.GT.MABS)GOTO 3 ASLEG=PLM2 GOTO 100 3 PLM1=(MABS+MABS+1)*X*PLM2 IF(L.GT.(MABS+1))GOTO 4 ASLEG=PLM1 GOTO 100 4 IU=MABS+2 DO 5 I=IU,L ASLEG=((I+I-1)*X*PLM1-(I+MABS-1)*PLM2)/(I-MABS) PLM2=PLM1 PLM1=ASLEG 5 CONTINUE 100 IF(M.GE.0)RETURN FAK=1.D0 IO=MABS+MABS DO 110 I=1,IO 110 FAK=FAK*(L-MABS+I) ASLEG=(-1)**MABS*ASLEG/FAK RETURN END c**************END OF PROGRAM FOR A-STATE********************