subroutine BKMP_pes(r,v) implicit real*8 (a-h,o-z) dimension r(3), dv(3) call h3tot_HHH_BKMP(r,v,dv,0) v=v-(-0.1744957296844722)+0.000001/27.2116 return end cmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmc cc-------------------------------------------------------------- c implicit double precision (a-h,o-z) c dimension x(18),Eknown(6),r(3),dVtot(3) c dimension dVdr(3),dVdrn(3) c parameter( npts=6, dr=1.d-08 ) c data x/ 1.757d0, 1.757d0, 3.514d0, c . 1.d0 , 2.d0 , 3.d0 , c . 3.d0 , 4.d0 , 5.d0 , c . 2.d0 , 2.d0 , 2.d0 , c . 1.d0 , 1.d0 , 1.d0 , c . 0.7d0 , 0.7d0 , 0.7d0 / c data Eknown/ -0.1592929999, -0.0928127174, -0.0561796335, c . -0.0739012770, 0.2122503352, 0.8559610809 / c call idh3(h3id) c write(6,*) ' H3 surface identification number: ',h3id c write(6,*) c write(6,*) ' comparison with known energies: ' c ideriv = 1 c write(6,6000) c do ipt=1,npts c ishift = (ipt-1)*3 c r(1) = x(ishift+1) c r(2) = x(ishift+2) c r(3) = x(ishift+3) c call h3tot(r,Vtot,dVtot,ideriv) c diff = Vtot - Eknown(ipt) c write(6,6005) r(1),r(2),r(3),Vtot,Eknown(ipt),diff c end do cc cc check analytical versus numerical derivatives for one geometry: cc save energy and analytical derivatives from last call c write(6,*) c write(6,*) 'compare analytical and numerical derivatives: ' c r(1) = 0.95d0 c r(2) = 1.05d0 c r(3) = 1.85d0 c call h3tot(r,Vref,dVdr,ideriv) c do i=1,3 c r(i) = r(i) + dr c call h3tot(r,Vtot,dVtot,ideriv) c dVdrn(i) = (Vtot-Vref)/dr c r(i) = r(i) - dr c end do c write(6,6015) r c write(6,6020) dVdr c write(6,6025) dVdrn c 6000 format(' r1 r2 r3 ', c . ' Vtot Eknown diff ') c 6005 format(3(f8.4,1x),2f16.10,2x,g12.6) c 6015 format(' r = ',3(1x,f12.6)) c 6020 format('analytical dV/dr = ',3(1x,e12.6)) c 6025 format('numerical dV/dr = ',3(1x,e12.6)) c stop c end cC ************************************************************************ ************************************************************************ c H3 surface evaluation routines c ============================== c block data BLKDT_HHH_BKMP c---------------------------------------------------------------------- c Note: vector 'ipass' is just used to pass debugging print flags around c if ipr.gt.0 ... print some details of various contributions to Vtot c if ipr.gt.1 ... also print out some details about derivatives c where ipr is set from various elements of ipass in various subroutines. c---------------------------------------------------------------------- common/pass_HHH_BKMP/ipass(100) data ipass/100*0/ end c subroutine idh3_HHH_BKMP(h3id) c----------------------------------------------------------------- c h3 surface identification number is 706 c----------------------------------------------------------------- implicit double precision (a-h,o-z) h3id = 706.0 return end c subroutine h3tot_HHH_BKMP(r,Vtot,dVtot,id) c----------------------------------------------------------------- c Calculate total H3 potential from all of its parts c if (id.gt.0) also calculate the derivative c all distances are in bohrs and all energies are in hartrees c c For a discussion of this surface, see: c A.I.Boothroyd, W.J.Keogh, P.G.Martin, M.R.Peterson c Journal of Chemical Physics 95 pp4343-4359 (Sept15/91) c Note: the surface parameter values as published lead to an anomolously c deep (several millihartree) van der Waals well for a very compact c H2 molecule (i.e., r < 0.8 bohr). After that paper was submitted, c this problem was fixed and the corrected Cbend coefficients are c used in this version of the surface (version 706); this change c should have negligible effect, since such small H2 molecules are c unbound. (The old coefficients are still in "subroutine vbcb" c but have been commented out.) c c any QUESTIONS/PROBLEMS/COMMENTS concerning this programme can be c addressed to: wkeogh@alchemy.chem.utoronto.ca c (if necessary, via boothroy@cita.utoronto.ca boothroy@utordop.bitnet c or pgmartin@cita.utoronto.ca pgmartin@utordop.bitnet ) c c---------------------------------------------------------------------- c version of july27/91 ... surf706d.out Cbend values put in c---------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension r(3),dVtot(3) dimension dVlon(3),dVas(3),dVbnda(3),dVbndb(3), . dCal(3), dCas(3),dCbnda(3),dCbndb(3) dimension dB1a(3),dB1b(3),dT(3),T(3) dimension dB2(3),dB3(6),expass(4) c {debugging:} dimension vb(2,25),cb(2,25),vbp(2,25),cbp(2,25) common/pass_HHH_BKMP/ipass(100) c ipr = ipass(1) c c FIRST, CHECK THAT IT IS A VALID H3 GEOMETRY, WITHIN TOLERANCE "DERR": c this test corrected on nov23/93 c generally, rlo + rmed > rhi c sometimes, rlo + rmed = rhi (linear geometry) c but if rlo + rmed < rhi there's a problem c derr = 1.d-5 rhi = dmax1(r(1),r(2)) if(r(3).ge.rhi)then rmid = rhi rhi = r(3) else rmid = dmax1( r(3) , dmin1(r(1),r(2)) ) endif rlo = dmin1(r(1),r(2),r(3)) cx if(rlo+rmid.gt.rhi+derr)then {old condition} if(rlo+rmid+derr.lt.rhi)then write(6,*) ' rlongest = ',rhi+derr write(6,*) ' rmiddle,rshortest = ',rmid,rlo stop ' STOP -- H3TOT Error: Not a valid H3 geometry. ' end if if(rlo.lt.0.2d0)then write(6,*) ' shortest distance: ',rlo coldversion stop ' STOP -- H3TOT Error: shortest distance < 0.2 bohr. ' end if c c zero everything to avoid any 'funny' values: Vlon = 0.d0 Vas = 0.d0 Vbnda = 0.d0 Vbndb = 0.d0 Cal = 0.d0 Cas = 0.d0 cbnda = 0.d0 cbndb = 0.d0 do i=1,3 dVlon(i) = 0.d0 dVas(i) = 0.d0 dVbnda(i)= 0.d0 dVbndb(i)= 0.d0 dCal(i) = 0.d0 dCas(i) = 0.d0 dCbnda(i)= 0.d0 dCbndb(i)= 0.d0 end do c zero the vb and cb arrays (used only for numerical derivatives) do i=1,25 vb(1,i) = 0.d0 vb(2,i) = 0.d0 cb(1,i) = 0.d0 cb(2,i) = 0.d0 end do call h3lond_HHH_BKMP( r, Vlon, dVlon ) call vascal_HHH_BKMP( r, Vas, dVas ) c now do any corrections required for compact geometries: call compac_HHH_BKMP(r,icompc,T,dT) c oct.3/90 compact routines only called for compact geometries: if(icompc.ge.1)then call csym_HHH_BKMP ( r, cal, dCal ) call casym_HHH_BKMP ( r, Cas, dCas, 1 ,T,dT) end if call VbCb_HHH_BKMP(r,icompc,T,dT,ipr, . Vbnda,Vbndb,dVbnda,dVbndb, . Cbnda,Cbndb,dCbnda,dCbndb) c add up the various parts of the potential: Vtot = Vlon + Vas + Vbnda + Vbndb + Cal + Cas + Cbnda + Cbndb c add up the various parts of the derivative: dVtot(1) = dVlon(1) + dVas(1) + dVbnda(1) + dVbndb(1) . + dCal(1) + dCas(1) + dCbnda(1) + dCbndb(1) dVtot(2) = dVlon(2) + dVas(2) + dVbnda(2) + dVbndb(2) . + dCal(2) + dCas(2) + dCbnda(2) + dCbndb(2) dVtot(3) = dVlon(3) + dVas(3) + dVbnda(3) + dVbndb(3) . + dCal(3) + dCas(3) + dCbnda(3) + dCbndb(3) if(ipr.gt.0) then write(7,*) write(7,*) 'h3tot_HHH_BKMP enter ----------------------------- $---' write(7,7000) r,icompc,vtot, . vlon,vas,cal,cas,vbnda,vbndb,cbnda,cbndb if(ipr.gt.1)then write(7,7100) dVtot,dVlon,dVas,dVbnda,dVbndb, . dCal,dCas,dCbnda,dCbndb end if write(7,7400) (vb(1,i),i=1,5),(vb(2,i),i=1,5) write(7,7500) (cb(1,i),i=1,9),(cb(2,i),i=1,9) write(7,7999) Vbnda,Vbndb,Cbnda,Cbndb, . dVbnda,dvbndb,dCbnda,dCbndb write(7,*) 'exiting subr.h3tot_HHH_BKMP' write(7,*) 'h3tot_HHH_BKMP exit ------------------------------ $---' end if return 7400 format('Vba values: ',5(1x,g12.6),/,'Vbb values: ',5(1x,g12.6)) 7500 format('Cba values: ',3(1x,f16.8),/, . ' ',3(1x,f16.8),/, . ' ',3(1x,f16.8),/, . 'Cbb values: ',3(1x,f16.8),/, . ' ',3(1x,f16.8),/, . ' ',3(1x,f16.8)) 7000 format(5x,' r =',3(1x,f16.10),' icompac = ',i1,/, . 5x,'Vtot = ',f18.12,/, . 5x,'Vlon = ',f18.12,' Vas = ',g18.12,/, . 5x,'Cal = ',g18.12,' Cas = ',g18.12,/, . 5x,'Vbnda = ',g18.12,' Vbndb = ',g18.12,/, . 5x,'Cbnda = ',g18.12,' Cbndb = ',g18.12) 7100 format(' dVtot = ',3(1x,g18.12),/, . ' dVlon = ',3(1x,g18.12),/, . ' dVas = ',3(1x,g18.12),/, . ' dVbnda = ',3(1x,g18.12),/, . ' dVbndb = ',3(1x,g18.12),/, . ' dCal = ',3(1x,g18.12),/, . ' dCas = ',3(1x,g18.12),/, . ' dCbnda = ',3(1x,g18.12),/, . ' dCbndb = ',3(1x,g18.12)) 7750 format(3(1x,g16.10)) 7751 format(/,3(1x,g16.10)) 7999 format(' Vbnda Vbndb Cbnda Cbndb ',/, . 4(e12.6,2x),/, . 'dVbnda: ',3(e16.10,1x),/, . 'dVbndb: ',3(e16.10,1x),/, . 'dCbnda: ',3(e16.10,1x),/, . 'dCbndb: ',3(e16.10,1x)) end c subroutine triplet_HHH_BKMP(r,E3,id) C-------------------------------------------------------------------- c version of oct4/90 surface626 values C H2 triplet curve and derivatives: c calculates triplet potential and first derivative c uses truhlar horowitz equation with our extension c uses the Johnson correction at short distances (r < rr) c if r .ge. rr use modified t/h triplet equation c if r .le. rl use the Johnson correction c in between use the transition equation C-------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension E3(3),E1(3) common/pass_HHH_BKMP/ipass(100) parameter( rl=0.95d0, rr=1.15d0 ) parameter( z1=1.d0, z2=2.d0, z4=4.d0) c triplet and johnson values from fit621: parameter(a1=-0.0253496194,a2=-29.2302126444, . a3=-50.7225015503,a4= 2.0452676876 , . a5=-12.2408908509,a6= 1.6733157383 ) parameter( C1=-0.4170298146519658, . C2=-0.0746027774843370,C3= 0.4297899952237434 ) c surface parameters from fit601.out c parameter(a1=-0.66129429,a2=-1.99434198,a3=-2.37604328, c . a4= 2.08107802,a5=-0.0313032510,a6=3.76546699, c . C1=-0.4222590135447196,C2=-0.0731117796738824, c . C3= 0.4295918082189010 ) ipr = ipass(9) E3(3) = 0.d0 if(r.ge.rr )then c modified truhlar/horowitz triplet equation: exdr = dexp( -a4*r ) rsq = r*r ra6 = r**(-a6) E3(1) = a1* ( a2 + r + a3*rsq + a5*ra6 )*exdr c first derivative of triplet curve: ra61 = r**(-a6-z1) E3(2) = a1*exdr* . ( z1 -a2*a4 +(z2*a3-a4)*r -a3*a4*rsq -a5*a6*ra61 -a4*a5*ra6 ) end if if(r.lt.rr )then dr = r- rl call vh2opt_HHH_BKMP(r,e1,2) if( r.le.rl ) then c Johnson triplet equation: E3(1) = E1(1) + c2*dr + c3 E3(2) = E1(2) + c2 else c Transition equation: E3(1) = E1(1) + c1*dr*dr*dr + c2*dr + c3 E3(2) = E1(2) + 3.d0*c1*dr*dr + c2 end if end if if(ipr.gt.0) write(7,7100) e3 7100 format(' triplet_HHH_BKMP: E3 = ',3(1x,f12.8)) return end C c subroutine vH2opt_HHH_BKMP(r,E,ideriv) C----------------------------------------------------------------- c version of jul 09/90 ... super duper speedy version C self-contained version of schwenke's H2 potential C all distances in bohrs and all energies in Hartrees C (1st deriv added on May 2 1989; 2nd deriv on May 28 1989) C----------------------------------------------------------------- implicit double precision (a-h,o-z) dimension E(3) common/pass_HHH_BKMP/ipass(100) parameter(a0=0.03537359271649620, a1= 2.013977588700072 , . a2= -2.827452449964767 , a3= 2.713257715593500 , . a4= -2.792039234205731 , a5= 2.166542078766724 , . a6= -1.272679684173909 , a7= 0.5630423099212294 , . a8= -0.1879397372273814 , a9= 0.04719891893374140 , . a10= -0.008851622656489644 , a11= 0.001224998776243630 , . a12= -1.227820520228028d-04 , a13= 8.638783190083473d-06, . a14= -4.036967926499151d-07 , a15= 1.123286608335365d-08, . a16= -1.406619156782167d-10 ) parameter( r0=3.5284882d0, DD=0.160979391d0, . c6=6.499027d0, c8=124.3991d0, c10=3285.828d0) c Ediss = 0.174445d0 Ediss = 0.174445d0*0.d0 E(1) = -999.d0 E(2) = -999.d0 E(3) = -999.d0 r2 = r *r r3 = r2 *r r4 = r3 *r r5 = r4 *r r6 = r5 *r r7 = r6 *r r8 = r7 *r r9 = r8 *r r10 = r9 *r r11 = r10*r r12 = r11*r r13 = r12*r r14 = r13*r r15 = r14*r r02 = r0*r0 r04 = r02*r02 r06 = r04*r02 rr2 = r2 + r02 rr4 = r4 + r04 rr6 = r6 + r06 rr25 = rr2*rr2*rr2*rr2*rr2 c general term: a(i)*r(i-1), i=0,16 alphaR = a0/r + a1 . + a2 *r + a3 *r2 + a4 *r3 + a5 *r4 + a6 *r5 . + a7 *r6 + a8 *r7 + a9 *r8 + a10*r9 + a11*r10 . + a12*r11 + a13*r12 + a14*r13 + a15*r14 + a16*r15 exalph = dexp(alphaR) vsr = DD*(exalph-1.d0)*(exalph-1.d0) - DD vlr = -C6/rr6 -C8/(rr4*rr4) -C10/rr25 E(1)= vsr + vlr c C calculate first derivative if required: if(ideriv.ge.1)then r3 = r2*r r5 = r4*r rr26 = rr25*rr2 rr43 = rr4*rr4*rr4 c general term: (i-1)*a(i)*r**(i-2) , i=0,16 dalphR = -a0/r2 + a2 + 2.d0*a3*r . + 3.d0 *a4 *r2 + 4.d0 *a5 *r3 + 5.d0 *a6 *r4 . + 6.d0 *a7 *r5 + 7.d0 *a8 *r6 + 8.d0 *a9 *r7 . + 9.d0 *a10*r8 + 10.d0 *a11*r9 + 11.d0 *a12*r10 . + 12.d0 *a13*r11 + 13.d0 *a14*r12 + 14.d0 *a15*r13 . + 15.d0 *a16*r14 dvsr = 2.d0 *DD *(exalph-1.d0) *exalph *dalphR dvlr = 6.d0*C6*r5 / (rr6*rr6) . + 8.d0*C8*r3 / rr43 + 10.d0*C10*r / rr26 E(2) = dvsr + dvlr end if c C calculate second derivative if required: if(ideriv.ge.2)then r10 = r6*r4 rr27 = rr26*rr2 rr44 = rr43*rr4 rr62 = rr6*rr6 rr63 = rr62*rr6 c general term: (i-1)*(i-2)*a(i)*r**(i-3), i=0,16 ddalph = 2.d0*a0/r3 + 2.d0*a3 + 6.d0*a4 *r . + 12.d0*a5 *r2 + 20.d0*a6 *r3 + 30.d0*a7 *r4 . + 42.d0*a8 *r5 + 56.d0*a9 *r6 + 72.d0*a10*r7 . + 90.d0*a11*r8 +110.d0*a12*r9 +132.d0*a13*r10 . +156.d0*a14*r11 +182.d0*a15*r12 +210.d0*a16*r13 ddvsr = 2.d0 *DD *exalph . *( (2.d0*exalph-1.d0)*dalphR*dalphR + (exalph-1.d0)*ddalph ) ddvlr =- 72.d0* C6 *r10 / rr63 -96.d0 *C8 *r6 / rr44 . -120.d0*C10 *r2 / rr27 +30.d0 *C6 *r4 / rr62 . + 24.d0* C8 *r2 / rr43 +10.d0 *C10 / rr26 E(3) = ddvsr + ddvlr end if return end c subroutine h3lond_HHH_BKMP(r,Vlon,dVlon) c---------------------------------------------------------------------- c version of may 12/90 ... derivatives corrected (0.5 changed to 0.25) c calculates the h3 london terms and derivatives c modified oct 7/89 to include eps**2 term which rounds off the c cusp in the h3 potential which occurs at equilateral triangle c configurations c---------------------------------------------------------------------- implicit double precision (a-h,o-z) real*8 Q(3),J(3),Jt dimension r(3),e1(3),e3(3),esing(3),etrip(3),dVlon(3) dimension de1(3),de3(3),dJ1(3),dJ2(3),dJ3(3) common/pass_HHH_BKMP/ipass(100) parameter( half=0.5d0, two=2.d0 , eps2=1.d-12 ) ipr=ipass(2) do i=1,3 call vh2opt_HHH_BKMP(r(i),esing,2) e1(i) = esing(1) de1(i) = esing(2) call triplet_HHH_BKMP(r(i),etrip,2) if(ipr.gt.0) write(7,7400) i,esing,etrip e3(i) = etrip(1) de3(i) = etrip(2) q(i) = half*(E1(i) + E3(i)) j(i) = half*(E1(i) - E3(i)) end do sumQ = q(1) + q(2) + q(3) sumJ = dabs( j(2)-j(1) )**2 . + dabs( j(3)-j(2) )**2 . + dabs( j(3)-j(1) )**2 Jt = half*sumJ + eps2 rootJt = dsqrt(Jt) Vlon = sumQ - rootJt if(ipr.gt.0) then write(7,7410) sumq,sumj write(7,7420) vlon,rootJt end if c calculate the derivatives with respect to r(i): dVlon(1) = half*(dE1(1)+de3(1)) . - 0.25d0*(two*j(1)-j(2)-j(3))*(dE1(1)-de3(1))/rootJt dVlon(2) = half*(dE1(2)+de3(2)) . - 0.25d0*(two*j(2)-j(3)-j(1))*(dE1(2)-de3(2))/rootJt dVlon(3) = half*(dE1(3)+de3(3)) . - 0.25d0*(two*j(3)-j(1)-j(2))*(dE1(3)-de3(3))/rootJt if(ipr.gt.0) then write(7,7000) r,e1,e3,vlon write(7,7100) q,j write(7,7200) dVlon end if 7000 format(' r = ',3(1x,f12.6),/, . ' E1 = ',3(1x,f12.8),/, . ' E3 = ',3(1x,f12.8),/, . ' vlon = ',1x,f12.8) 7100 format(13x,'q = ',3(1x,e12.6),/,13x,'j = ',3(1x,e12.6)) 7200 format(' dVlon = ',3(1x,g12.6)) 7400 format('from subr.london: ',/, . ' using r',i1,': Esinglet=',3(1x,f12.8),/, . ' ',1x,' Etriplet=',3(1x,f12.8)) 7410 format(' sumQ = ',g12.6,' sumJ = ',g12.6) 7420 format(' Vlon = ',f12.8,' rootJt = ',g12.6) return end c subroutine Vascal_HHH_BKMP(rpass,Vas,dVas) c------------------------------------------------------------------ c version of oct11/90 fit632c.out values c version of oct5/90 surf626 values c calculate the asymmetric correction term and its derivatives c see equations [14] to [16] of truhlar/horowitz 1978 paper c------------------------------------------------------------------ implicit double precision (a-h,o-z) common/pass_HHH_BKMP/ipass(100) dimension dVas(3),dA(3),dS(3),rpass(3) c--- Vasym values from fit632c ---------------------- parameter( . aa1=0.3438222224E-02,aa2=0.1398145763E-02, . aa3=-.1923999449E-03,aa4=0.9712737075E-05, . aa5=-.1263794562E-06,aa6=0.5181432712E+00, . aa7=-.9487002995E-03 ) ipr = ipass(3) r1 = rpass(1) r2 = rpass(2) r3 = rpass(3) R = r1 + r2 + r3 Rsq = R*R Rcu = Rsq*R C calculate the Vas term first (eq.14 of truhlar/horowitz) call acalc_HHH_BKMP(r1,r2,r3,A,dA) A2 = A *A A3 = A2*A A4 = A3*A A5 = A4*A exp1 = dexp(-aa1*Rcu) exp6 = dexp(-aa6*R) S = aa2*A2 + aa3*A3 + aa4*A4 + aa5*A5 Vas = S*exp1 + aa7*A2 *exp6 / R dS(1) = ( 2.d0*aa2*A +3.d0*aa3*A2 . +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(1) dVas(1) = -3.d0*aa1*Rsq*S*exp1 + dS(1)*exp1 . -aa7*A2*exp6/Rsq + 2.d0*aa7*A*dA(1)*exp6/R . -aa6*aa7*A2*exp6/R dS(2) = ( 2.d0*aa2*A +3.d0*aa3*A2 . +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(2) dVas(2) = -3.d0*aa1*Rsq*S*exp1 + dS(2)*exp1 . -aa7*A2*exp6/Rsq + 2.d0*aa7*A*dA(2)*exp6/R . -aa6*aa7*A2*exp6/R dS(3) = ( 2.d0*aa2*A +3.d0*aa3*A2 . +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(3) dVas(3) = -3.d0*aa1*Rsq*S*exp1 + dS(3)*exp1 . -aa7*A2*exp6/Rsq + 2.d0*aa7*A*dA(3)*exp6/R . -aa6*aa7*A2*exp6/R if(ipr.gt.0) write(7,7100) Vas,dS,dVas 7100 format(' from subr.vascal: Vas = ', 1x,g12.6, /, . ' dS = ',3(1x,g12.6),/, . ' dVas = ',3(1x,g12.6)) return end c subroutine acalc_HHH_BKMP(r1,r2,r3,A,dA) c--------------------------------------------------------------------- c version of may 1, 1990 c based on equations from notes date april 4, 1990 c calculate the A value and its derivatives c A = dabs[ (r1-r2)*(r2-r3)*(r3-r1) ] c--------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension dA(3) common/pass_HHH_BKMP/ipass(100) ipr = ipass(11) A = (r1-r2)*(r2-r3)*(r3-r1) dA(1) = ( -2.d0*r1 + r2 + r3 )*(r2-r3) dA(2) = ( -2.d0*r2 + r3 + r1 )*(r3-r1) dA(3) = ( -2.d0*r3 + r1 + r2 )*(r1-r2) if(A.lt.0.d0)then A = -A dA(1) = -dA(1) dA(2) = -dA(2) dA(3) = -dA(3) end if if(ipr.gt.0) write(7,7100) A,dA 7100 format('Acalc enter ---------------------------------------',/, . ' A = ',f15.10,' dA = ',3(1x,g12.6),/, . 'Acalc exit ----------------------------------------') return end c subroutine compac_HHH_BKMP(r,icompc,T,dT) C--------------------------------------------------------------- c version of may 14/90 ... calculates T and dT values also c decide whether or not this particular geometry is compact, c that is, are any of the three distances smaller than the c rr value from the johnson correction. c----------------------------------------------------------------- implicit double precision (a-h,o-z) dimension r(3),T(3),dT(3) common/pass_HHH_BKMP/ipass(100) parameter( rr = 1.15d0, rp = 1.25d0 ) c first see if this is a compact geometry: do i=1,3 T(i) = 0.d0 dT(i)= 0.d0 end do icompc = 0 if(r(1).lt.rr) icompc = icompc + 1 if(r(2).lt.rr) icompc = icompc + 1 if(r(3).lt.rr) icompc = icompc + 1 if(icompc.eq.0) return c calculate the T(i) values: do i=1,3 if(r(i).lt.rr) then top = rr-r(i) bot = rp-r(i) top2 = top * top top3 = top2 * top bot2 = bot * bot T(i) = top3/bot dT(i) = -3.d0*top2/bot + top3/bot2 end if end do return end C subroutine casym_HHH_BKMP ( r,Cas,dCas,id,T,dT ) c--------------------------------------------------------------- c version of sept14/90 c the compact asymmetric correction term and derivatives c--------------------------------------------------------------- implicit double precision (a-h,o-z) dimension r(3),T(3) dimension dPr(3),dT(3),dsumT(3),dterm1(3),determ(3),dseries(3), . dCas(3),dA(3) common/pass_HHH_BKMP/ipass(100) c---- Casym values from fit633a parameter( . u1=0.1537481166E+00, u2=0.2745950036E+00, u3=0.7501206780E-02, . u4=0.3119136023E+01, u5=0.9969170798E+00, u6=-.3373682823E+01, . u7=0.6807215913E+00, u8=-.4920325491E-01, u9=-.3919467989E+01, .u10=0.5085532326E+01,u11=-.1415264778E+01,u12=0.1138681785E+00, .u13=-.1525367566E-02) ipr = ipass(6) cas = 0.d0 dCas(1) = 0.d0 dCas(2) = 0.d0 dCas(3) = 0.d0 call acalc_HHH_BKMP(r(1),r(2),r(3),A,dA) A2 = A*A c if(A.eq.0.d0) return sumT = T(1) + T(2) + T(3) Sr = r(1) + r(2) + r(3) Pr = r(1) * r(2) * r(3) Sr2 = Sr*Sr Sr3 = Sr2*Sr Pr2 = Pr*Pr Pr3 = Pr2*Pr c write out the series explicitly: series = 1.d0 + u4/Pr2 + u5/Pr + u6 + u7*Pr + u8*Pr2 . + A*(u9/Pr2 + u10/Pr + u11 + u12*Pr + u13*Pr2) term1 = u1/Pr**u2 eterm = dexp(-u3*Sr3) Cas = sumt * a2 * term1 * series * eterm if(ipr.gt.0) then write(7,*) write(7,*) 'Casym: ------------------' write(7,*) ' T(i) = ',T write(7,*) 'dT(i) = ',dT write(7,*) ' id = ',id write(7,7400) series,term1,eterm,Cas write(7,*) ' Cas = ',Cas end if c {calculate the derivative as well:} if(id.gt.0)then dPr(1) = r(2)*r(3) dPr(2) = r(3)*r(1) dPr(3) = r(1)*r(2) do i=1,3 dsumt(i) = dt(i) dterm1(i)= -1.d0*u1*u2*Pr**(-u2-1.d0)*dPr(i) determ(i)= eterm *(-3.d0*u3*Sr2) dseries(i)= . dPr(i)*( -2.d0*u4/Pr3 -u5/Pr2 +u7 +2.d0*u8*Pr ) . +dA(i)*( u9/Pr2 +u10/Pr +u11 +u12*Pr +u13*Pr2 ) . +A*dPr(i)*( -2.d0*u9/Pr3 -u10/Pr2 +u12 + 2.d0*u13*Pr ) dCas(i) = . dsumt(i) * a2 * term1 * series * eterm . + 2.d0*a*dA(i) * sumt * term1 * series * eterm . + dterm1(i) * sumt * a2 * series * eterm . + dseries(i) * sumt * a2 * term1 * eterm . + determ(i) * sumt * a2 * term1 * series end do if(ipr.gt.0) then write(7,7500) dPr,dt,dterm1,determ,dseries,dCas write(7,*) ' dCas = ',dCas end if end if return 7400 format('subr.Cas: series = ',g12.6,' term1 = ',g12.6,/, . ' eterm = ',g12.6,' Cas = ',g12.6) 7500 format(' derivatives: dPr = ',3(1x,g12.6),/, . ' dT = ',3(1x,g12.6),/, . ' dterm1 = ',3(1x,g12.6),/, . ' determ = ',3(1x,g12.6),/, . ' dseries = ',3(1x,g12.6),/, . ' dCas = ',3(1x,g12.6)) end c subroutine csym_HHH_BKMP(r,Cal,dCal) c--------------------------------------------------------------- c version of oct12/90 surf636 values c calculate the 'compact all' correction term and derivatives c a correction term (added sept 11/89), which adds a small c correction to all compact geometries c--------------------------------------------------------------- implicit double precision (a-h,o-z) dimension r(3),t(3),dCal(3),G(3) dimension dG(3),sumv(3) common/pass_HHH_BKMP/ipass(100) parameter( rr=1.15d0, rp=1.25d0 ) c---- Csym values from fit633a parameter( . v1=-.2070776049E+00, v2=-.5350898737E+00, v3=0.1011861942E-01) cal = 0.d0 ipr = ipass(5) Sr = r(1)+r(2)+r(3) Sr2 = Sr*Sr Sr3 = Sr*Sr2 exp3 = dexp( -v3*Sr3 ) dexp3 = -3.d0*v3*Sr2*exp3 do i=1,3 ri = r(i) rrri = rr-ri rrri2 = rrri*rrri rrri3 = rrri*rrri2 rpri = rp-ri rpri2 = rpri*rpri G(i) = 0.d0 dG(i) = 0.d0 sumv(i) = v1+v1*v2*ri if(ri.lt.rr) then G(i) = (rrri3/rpri) *sumv(i) dG(i) = (rrri3/rpri2)*sumv(i) . -3.d0*(rrri2/rpri)*sumv(i) . + (rrri3/rpri)*v1*v2 end if end do sumG = G(1) + G(2) + G(3) cal = sumG*exp3 dCal(1) = dG(1)*exp3 + sumG*dexp3 dCal(2) = dG(2)*exp3 + sumG*dexp3 dCal(3) = dG(3)*exp3 + sumG*dexp3 if(ipr.gt.0)then write(7,*) write(7,7100) Cal write(7,7000) G,dG,sumv,dCal write(7,*) 'exiting subr.csym_HHH_BKMP' end if 7100 format('Csym enter ----------------------------------------',/, . ' Cal = ',1x,g20.14) 7000 format(' G = ',3(1x,g12.6),/, . ' dG = ',3(1x,g12.6),/, . ' sumv = ',3(1x,g12.6),/, . ' dCal = ',3(1x,g16.10),/, . 'Csym exit -----------------------------------------') return end c subroutine VbCb_HHH_BKMP(rpass,icompc,T,dT,ipr, . Vbnda,Vbndb,dVbnda,dVbndb, . Cbnda,Cbndb,dCbnda,dCbndb) C--------------------------------------------------------------- c jul24/91 fit705 cbend values added c feb27/91 modified to match equation in h3 paper more closely c nov 4/90 Vbend coefficients now a,g Cbend coeff's still c,d c in this version, the derivatives are always calculated c subroutines vbend, cbend and B1ab all combined into this one c module in order to improve efficiency by not having to pass c around B1a,B1b,B2,B3a,B3b functions and derivatives c B1a = 1 - sum of [ P1(cos(theta(i))) ] c B1b = 1 - sum of [ P3(cos(theta(i))) ] c-------------------------------------------------------------- implicit double precision (a-h,o-z) c arrays required for B1,B2,B3 calculations: dimension rpass(3),dB1a(3),dB1b(3),th(3), . dB2(3),dB3(6) c arrays required for Vbend calculations: dimension dVbnda(3),dVbndb(3) dimension dVb1(3),dvb2(3),dvb3(3),dvb4(3),dvb5(3),dVbnd(3) dimension Vbnd(2) c {passed out for debugging only:} dimension vb(2,25) c arrays required for Cbend calculations: dimension Cbnds(2),dCbnds(2,3) dimension dCb1(3),dCb2(3),dCb3(3),dCb4(3),dCb5(3),dCb6(3), . dCb7(3),dCb8(3) dimension T(3),dP(3) dimension dT(3),dCbnda(3),dCbndb(3) c {passed out for debugging only:} dimension cb(2,25) c {3/8 and 5/8:} parameter( z58=0.625d0, z38=0.375d0 ) c parameters required for Vbend calculations: parameter(beta1=0.52d0, beta2=0.052d0, beta3=0.79d0 ) c-- Vbenda values from fit634b parameter( .a11=-.2918252280E+03,a12=0.2164569141E+02,a13=-.4005497543E+00, .a21=-.2890774947E+01,a22=0.1032032542E+02,a23=0.2681748056E+02, .a24=0.2633751055E+02,a31=0.6180641351E+01,a32=0.5037667041E+01, .a41=-.1125570079E+00,a42=-.3176529304E-02,a43=0.9068915355E+00, .a44=-.7228418516E+00,a51=0.2785898232E+03,a52=0.4764442446E+02, .a53=0.8621545662E+01) c-- Vbendb values from fit634b parameter( .g11=-.4241912309E+02,g12=0.2951365281E+01,g13=-.4840201514E-01, .g21=-.1159168549E-03,g22=0.8688003567E-02,g23=-.3486923900E-02, .g24=0.8312212839E-03,g31=0.5621810473E-01,g32=-.9776930747E-02, .g41=-.1178456251E-01,g42=0.3491086729E-02,g43=0.7430516993E-01, .g44=-.9643636957E-01,g51=0.4735533782E+02,g52=0.3001038808E+01, .g53=0.1896630453E+01) c parameters required for Cbend calculations: c-- Cbenda values from fit702c ------------------- cx parameter( cx .c11=0.3408013485E+04,c12=-.2569675505E+04,c13=-.4171156109E+03, cx .c14=-.6179297444E+04,c15=0.3235223485E+04,c21=-.3409049347E+03, cx .c22=0.9612365597E+02,c23=0.4104361425E+03,c24=0.6238298637E+03, cx .c31=0.2655770314E+04,c32=0.5314085916E+03,c41=0.2222002477E+02, cx .c42=-.6611541362E+01,c43=-.2541255632E+03,c44=0.7752854954E+02, cx .c51=0.9613217362E+04,c52=-.7010043704E+03,c53=0.1406657509E+03, cx .c61=-.1677846467E+04,c62=0.7411196008E+04,c63=-.1104074171E+05, cx .c71=-.1077399844E+03,c72=-.6598646429E+01,c73=0.2035775472E+00, cx .c74=0.1097670442E+04,c75=0.4873456879E+02,c81=0.1005135589E+03, cx .c82=-.1972894789E+02,c83=-.8495015577E+04,c84=-.6106570657E+02) c-- Cbendb values from fit702c ------------------- cx parameter( cx .d11=0.8450990097E+03,d12=0.1466587231E+03,d13=0.1239413655E+02, cx .d14=-.5649327722E+03,d15=0.3112693920E+03,d21=0.8029625257E-01, cx .d22=0.2379962852E-01,d23=0.1233566879E+00,d24=-.3582414079E-01, cx .d31=-.1012436386E+03,d32=0.8904152887E-02,d41=0.8740633605E+01, cx .d42=-.2937137124E+01,d43=-.5497174488E+02,d44=0.1656850593E+02, cx .d51=0.3859496208E+05,d52=0.8458336162E+02,d53=-.1431690434E+02, cx .d61=0.5487645288E+02,d62=-.1954058540E+03,d63=-.7226235345E+02, cx .d71=-.1895543609E+02,d72=0.1281231795E+00,d73=0.1291641223E+00, cx .d74=-.3052282568E+03,d75=-.6194562183E+01,d81=-.1746535079E+02, cx .d82=-.8402855300E+00,d83=-.3878264141E+05,d84=0.5655392575E+01) c-- Cbenda values from fit705a.out ------------- cxx parameter( cxx .c11=0.7138332997E+04,c12=-.3941467854E+04,c13=0.1214839191E+03, cxx .c14=-.9715171923E+04,c15=0.4549112075E+04,c21=-.7153573233E+02, cxx .c22=0.1410801562E+03,c23=0.5310286989E+03,c24=0.6632669305E+03, cxx .c31=0.9648491179E+03,c32=0.1336994724E+03,c41=0.3986106227E+01, cxx .c42=-.2071394753E+01,c43=-.1304296435E+03,c44=0.4433517239E+02, cxx .c51=-.2160974730E+05,c52=-.5040404494E+03,c53=0.1150251071E+03, cxx .c61=-.1900650362E+04,c62=0.8065979973E+04,c63=-.1022333250E+05, cxx .c71=-.4023633922E+02,c72=-.5593869000E+01,c73=0.2901552041E+00, cxx .c74=0.1332299065E+04,c75=0.4830933142E+02,c81=0.1527546398E+03, cxx .c82=-.2277696550E+02,c83=0.2214880233E+05,c84=-.5951714583E+02) c-- Cbendb values from fit705a.out ------------- cxx parameter( cxx .d11=0.1733880674E+04,d12=-.5105474641E+03,d13=0.5739180116E+02, cxx .d14=-.1549259263E+04,d15=0.6494865366E+03,d21=-.1180847496E+00, cxx .d22=0.3847529877E-01,d23=0.1035799123E+00,d24=-.2783059123E-01, cxx .d31=-.4049790116E+02,d32=0.3623102614E-01,d41=0.2132254354E+01, cxx .d42=-.1337841873E+01,d43=-.2911759497E+02,d44=0.9698950726E+01, cxx .d51=0.7851192036E+05,d52=0.1244750132E+03,d53=-.2157749426E+02, cxx .d61=0.1945719436E+02,d62=-.4336295493E+02,d63=0.8025442515E+02, cxx .d71=0.5232787178E+01,d72=-.2082854432E+01,d73=0.1740967673E+00, cxx .d74=-.3069559809E+02,d75=0.9099007208E+01,d81=-.6088433744E+01, cxx .d82=-.9840250626E+00,d83=-.7880550115E+05,d84=-.9084945347E+01) c-- Cbenda values from fit706d.out parameter( .c11=0.7107064647E+04,c12=-.3728421267E+04,c13=0.1757654624E+03, .c14=-.9725998132E+04,c15=0.4665086074E+04,c21=-.4435165986E+02, .c22=0.1604477309E+03,c23=0.5805142925E+03,c24=0.6892349445E+03, .c31=0.6581730442E+03,c32=0.6078389739E+02,c41=0.2885182566E+01, .c42=-.1728169916E+01,c43=-.1119535503E+03,c44=0.4052536250E+02, .c51=0.2540505673E+03,c52=-.5762083627E+03,c53=0.1295901320E+03, .c61=-.2131706075E+04,c62=0.9084452020E+04,c63=-.1138253963E+05, .c71=-.3964298833E+02,c72=-.5019979693E+01,c73=0.2906541488E+00, .c74=0.1212943686E+04,c75=0.4140463415E+02,c81=0.1752855549E+03, .c82=-.2496320107E+02,c83=0.3765413052E+03,c84=-.5480488130E+02) c-- Cbendb values from fit706d.out parameter( .d11=0.1917166552E+04,d12=-.6542563392E+03,d13=0.6793758367E+02, .d14=-.1694968847E+04,d15=0.6866649703E+03,d21=-.2137567948E+00, .d22=0.4975938228E-01,d23=0.9364998295E-01,d24=-.2444320779E-01, .d31=-.2863126914E+02,d32=0.5443219625E-01,d41=0.9673956120E+00, .d42=-.1160159706E+01,d43=-.2424199759E+02,d44=0.8569424490E+01, .d51=-.6517635862E+04,d52=0.1518098147E+03,d53=-.2706514366E+02, .d61=0.4308392956E+02,d62=-.1234851732E+03,d63=0.2320626055E+03, .d71=0.1049541418E+02,d72=-.2424169341E+01,d73=0.1745646946E+00, .d74=0.2603615561E+02,d75=0.1345799970E+02,d81=-.6653710513E+01, .d82=-.2576854447E+00,d83=0.6172608425E+04,d84=-.1328142473E+02) c feb27/91 new c51 = c51+c83; new d51 = d51+d83 cx1 = c51 + c83 dx1 = d51 + d83 r1 = rpass(1) r2 = rpass(2) r3 = rpass(3) t1 = r1*r1 -r2*r2 -r3*r3 t2 = r2*r2 -r3*r3 -r1*r1 t3 = r3*r3 -r1*r1 -r2*r2 c calculate the cosines of the three internal angles: c1 = t1 / (-2.d0*r2*r3) c2 = t2 / (-2.d0*r3*r1) c3 = t3 / (-2.d0*r1*r2) sum = c1 + c2 + c3 B1a = 1.d0 - sum c1cube = c1*c1*c1 c2cube = c2*c2*c2 c3cube = c3*c3*c3 cos3t1 = 4.d0*c1cube - 3.d0*c1 cos3t2 = 4.d0*c2cube - 3.d0*c2 cos3t3 = 4.d0*c3cube - 3.d0*c3 sumb = cos3t1 + cos3t2 + cos3t3 B1b = 1.d0 - ( z58*sumb + z38*sum ) c calculate derivatives if desired dc1dr1 = -r1/(r2*r3) dc2dr2 = -r2/(r1*r3) dc3dr3 = -r3/(r1*r2) dc1dr2 = ( t1/(r2*r2) + 2.d0 )/(2.d0*r3) dc1dr3 = ( t1/(r3*r3) + 2.d0 )/(2.d0*r2) dc2dr1 = ( t2/(r1*r1) + 2.d0 )/(2.d0*r3) dc2dr3 = ( t2/(r3*r3) + 2.d0 )/(2.d0*r1) dc3dr1 = ( t3/(r1*r1) + 2.d0 )/(2.d0*r2) dc3dr2 = ( t3/(r2*r2) + 2.d0 )/(2.d0*r1) db1a(1) = -1.d0*( dc1dr1 + dc2dr1 + dc3dr1 ) db1a(2) = -1.d0*( dc1dr2 + dc2dr2 + dc3dr2 ) db1a(3) = -1.d0*( dc1dr3 + dc2dr3 + dc3dr3 ) d1 = 12.d0*c1*c1 - 3.d0 d2 = 12.d0*c2*c2 - 3.d0 d3 = 12.d0*c3*c3 - 3.d0 db1b(1) = -z58*( d1*dc1dr1 + d2*dc2dr1 + d3*dc3dr1 ) . -z38*( dc1dr1 + dc2dr1 + dc3dr1 ) db1b(2) = -z58*( d1*dc1dr2 + d2*dc2dr2 + d3*dc3dr2 ) . -z38*( dc1dr2 + dc2dr2 + dc3dr2 ) db1b(3) = -z58*( d1*dc1dr3 + d2*dc2dr3 + d3*dc3dr3 ) . -z38*( dc1dr3 + dc2dr3 + dc3dr3 ) d if(ipr.gt.1)then d write(7,*) d write(7,*) ' from subr.vbcb -----------------------------' d write(7,6000) rpass d write(7,6010) b1a,b1b d write(7,6020) db1a,db1b d if(ipr.gt.2)then d write(7,7000) dc1dr1,dc1dr2,dc1dr3, d . dc2dr1,dc2dr2,dc2dr3, d . dc3dr1,dc3dr2,dc3dr3 d write(7,7010) t1,t2,t3,c1,c2,c3,d1,d2,d3 d end if d end if c c calculate the quantities used by both Vbend and Cbend: c R = r1 + r2 + r3 Rsq = R*R B2 = 1.d0/r1 + 1.d0/r2 + 1.d0/r3 B3 = (r2-r1)*(r2-r1) + (r3-r2)*(r3-r2) + (r1-r3)*(r1-r3) eps2 = 1.d-12 B3b = dsqrt( B3 + eps2 ) dB2(1) = -1.d0/(r1*r1) dB2(2) = -1.d0/(r2*r2) dB2(3) = -1.d0/(r3*r3) dB3(1) = 4.d0*r1 -2.d0*r2 -2.d0*r3 dB3(2) = 4.d0*r2 -2.d0*r3 -2.d0*r1 dB3(3) = 4.d0*r3 -2.d0*r1 -2.d0*r2 dB3(4) = 0.5d0*dB3(1)/B3b dB3(5) = 0.5d0*dB3(2)/B3b dB3(6) = 0.5d0*dB3(3)/B3b exp1 = dexp(-beta1*R) exp2 = dexp(-beta2*Rsq) exp7 = dexp(-beta3*R) dexp1 = -beta1*exp1 dexp2 = -2.d0*beta2*R*exp2 dexp7 = -beta3*exp7 c do the vbnda calculations: B1 = B1a B12 = B1*B1 B13 = B12*B1 B14 = B13*B1 B15 = B14*B1 asum = a11 + a12*R + a13*Rsq bsum = a21*B12 + a22*B13 + a23*B14 + a24*B15 csum = a31*B1*exp1 + a32*B12*exp2 dsum1 = a41*exp1 + a42*exp2 dsum2 = a43*exp1 + a44*exp2 fsum = a51 + a52*R + a53*Rsq vb(1,1) = B1*asum*exp1 vb(1,2) = bsum * exp2 vb(1,3) = B2*csum vb(1,4) = B1*B3*dsum1 + B1*B3b*dsum2 vb(1,5) = B1* fsum *exp7 Vbnd(1) = vb(1,1)+vb(1,2)+vb(1,3)+vb(1,4)+vb(1,5) c dasum = a12+2.d0*a13*R dbsum = 2.d0* a21*B1 + 3.d0* a22*B12 . + 4.d0* a23*B13 + 5.d0* a24*B14 ddsum1 = a41*dexp1 + a42*dexp2 ddsum2 = a43*dexp1 + a44*dexp2 dfsum = a52 +2.d0*a53*R dVb1(1) = dB1a(1)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1 dVb1(2) = dB1a(2)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1 dVb1(3) = dB1a(3)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1 c dVb2(1) = dbsum*dB1a(1)*exp2 + bsum*dexp2 dVb2(2) = dbsum*dB1a(2)*exp2 + bsum*dexp2 dVb2(3) = dbsum*dB1a(3)*exp2 + bsum*dexp2 c c calculate the Vb3 derivatives: dVb3(1)=dB2(1)*csum + B2*( a31*dB1a(1)*exp1 + a31*B1*dexp1 . +2.d0* a32*B1*dB1a(1)*exp2 + a32*B12*dexp2) dVb3(2)=dB2(2)*csum + B2*( a31*dB1a(2)*exp1 + a31*B1*dexp1 . +2.d0* a32*B1*dB1a(2)*exp2 + a32*B12*dexp2) dVb3(3)=dB2(3)*csum + B2*( a31*dB1a(3)*exp1 + a31*B1*dexp1 . +2.d0* a32*B1*dB1a(3)*exp2 + a32*B12*dexp2) c c calculate the Vb4 derivatives (may 27/90): dVb4(1)= dB1a(1)*B3 *dsum1+ B1*dB3(1)*dsum1 + B1*B3 *ddsum1 . + dB1a(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2 dVb4(2)= dB1a(2)*B3 *dsum1+ B1*dB3(2)*dsum1 + B1*B3 *ddsum1 . + dB1a(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2 dVb4(3)= dB1a(3)*B3 *dsum1+ B1*dB3(3)*dsum1 + B1*B3 *ddsum1 . + dB1a(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2 c c calculate the Vb5 derivatives: dVb5(1) = dB1a(1)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7 dVb5(2) = dB1a(2)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7 dVb5(3) = dB1a(3)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7 c c calculate the overall derivatives: dVbnda(1) = dvb1(1)+dvb2(1)+dvb3(1)+dvb4(1)+dvb5(1) dVbnda(2) = dvb1(2)+dvb2(2)+dvb3(2)+dvb4(2)+dvb5(2) dVbnda(3) = dvb1(3)+dvb2(3)+dvb3(3)+dvb4(3)+dvb5(3) Vbnda = Vbnd(1) c c---- now do the vbendb calculations -------- B1 = B1b B12 = B1*B1 B13 = B12*B1 B14 = B13*B1 B15 = B14*B1 asum = g11 + g12*R + g13*Rsq bsum = g21*B12 + g22*B13 + g23*B14 + g24*B15 csum = g31*B1*exp1 + g32*B12*exp2 dsum1 = g41*exp1 + g42*exp2 dsum2 = g43*exp1 + g44*exp2 fsum = g51 + g52*R + g53*Rsq vb(2,1) = B1*asum*exp1 vb(2,2) = bsum * exp2 vb(2,3) = B2*csum vb(2,4) = B1*B3*dsum1 + B1*B3b*dsum2 vb(2,5) = B1* fsum *exp7 Vbnd(2) = vb(2,1)+vb(2,2)+vb(2,3)+vb(2,4)+vb(2,5) c dasum = g12+2.d0*g13*R dbsum = 2.d0* g21*B1 + 3.d0* g22*B12 . +4.d0* g23*B13 + 5.d0* g24*B14 ddsum1 = g41*dexp1 + g42*dexp2 ddsum2 = g43*dexp1 + g44*dexp2 dfsum = g52 +2.d0*g53*R dVb1(1) = dB1b(1)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1 dVb1(2) = dB1b(2)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1 dVb1(3) = dB1b(3)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1 c dVb2(1) = dbsum*dB1b(1)*exp2 + bsum*dexp2 dVb2(2) = dbsum*dB1b(2)*exp2 + bsum*dexp2 dVb2(3) = dbsum*dB1b(3)*exp2 + bsum*dexp2 c c calculate the Vb3 derivatives: dVb3(1)=dB2(1)*csum + B2*( g31*dB1b(1)*exp1 + g31*B1*dexp1 . +2.d0* g32*B1*dB1b(1)*exp2 + g32*B12*dexp2) dVb3(2)=dB2(2)*csum + B2*( g31*dB1b(2)*exp1 + g31*B1*dexp1 . +2.d0* g32*B1*dB1b(2)*exp2 + g32*B12*dexp2) dVb3(3)=dB2(3)*csum + B2*( g31*dB1b(3)*exp1 + g31*B1*dexp1 . +2.d0* g32*B1*dB1b(3)*exp2 + g32*B12*dexp2) c c calculate the Vb4 derivatives (may 27/90): dVb4(1)= dB1b(1)*B3 *dsum1+ B1*dB3(1)*dsum1 + B1*B3 *ddsum1 . + dB1b(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2 dVb4(2)= dB1b(2)*B3 *dsum1+ B1*dB3(2)*dsum1 + B1*B3 *ddsum1 . + dB1b(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2 dVb4(3)= dB1b(3)*B3 *dsum1+ B1*dB3(3)*dsum1 + B1*B3 *ddsum1 . + dB1b(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2 c c calculate the Vb5 derivatives: dVb5(1) = dB1b(1)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7 dVb5(2) = dB1b(2)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7 dVb5(3) = dB1b(3)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7 c c calculate the overall derivatives: dVbndb(1) = dvb1(1)+dvb2(1)+dvb3(1)+dvb4(1)+dvb5(1) dVbndb(2) = dvb1(2)+dvb2(2)+dvb3(2)+dvb4(2)+dvb5(2) dVbndb(3) = dvb1(3)+dvb2(3)+dvb3(3)+dvb4(3)+dvb5(3) Vbndb = Vbnd(2) c c now calculate the Cbend correction terms: c Cbnda uses the B1a formula and Cbndb uses the B1b formula if(icompc.eq.0) return sumT = T(1) + T(2) + T(3) Rcu = Rsq*R P = r1*r2*r3 Psq = P*P Pcu = Psq*P dP(1) = r2*r3 dP(2) = r3*r1 dP(3) = r1*r2 Cbnds(1) = 0.d0 Cbnds(2) = 0.d0 dCbnda(1) = 0.d0 dCbnda(2) = 0.d0 dCbnda(3) = 0.d0 dCbndb(1) = 0.d0 dCbndb(2) = 0.d0 dCbndb(3) = 0.d0 c calculate exponentials and derivatives (common to cbnda and cbndb): exp7 = dexp( -beta2*Rcu ) dexp7 = -3.d0*beta2*Rsq*exp7 c ----- calculate the Cbnda correction term ---------------- B1 = B1a B12 = B1*B1 B13 = B12*B1 B14 = B13*B1 B15 = B14*B1 asum = c11 + c12*R + c13*Rsq + c14/R + c15/Rsq bsum = c21*B12 + c22*B13 + c23*B14 + c24*B15 csum = c31*B1 *exp1 + c32*B12*exp2 dsum1 = c41*exp1 + c42*exp2 dsum2 = c43*exp1 + c44*exp2 fsum = cx1 + c52*R + c53*Rsq gsum = c61 + c62/R + c63/Rsq aasum = c71 + c72*P + c73*Psq + c74/P + c75/Psq ffsum = c81*P + c82*Psq + c84/Psq dasum = c12 + 2.d0*c13*R -c14/Rsq -2.d0*c15/Rcu dbsum = 2.d0*c21*B1 +3.d0*c22*B12 +4.d0*c23*B13 +5.d0*c24*B14 ddsum1= c41*dexp1 + c42*dexp2 ddsum2= c43*dexp1 + c44*dexp2 dfsum = c52 + 2.d0*c53*R dgsum = -c62/Rsq - 2.d0*c63/Rcu daasum = c72 +2.d0*c73*P -c74/Psq -2.d0*c75/Pcu dffsum = c81 + 2.d0*c82*P -2.d0*c84/Pcu cb(1,1) = B1*asum * exp1 /P cb(1,2) = bsum * exp2 cb(1,3) = B2*csum cb(1,4) = B1*B3*dsum1 + B1*B3b*dsum2 cb(1,5) = B1 * fsum * exp7 /P cb(1,6) = B1 * gsum * exp7 cb(1,7) = B1 * aasum * exp2 cb(1,8) = B1 * ffsum * exp7 Cbnds(1) = cb(1,1)+cb(1,2)+cb(1,3)+cb(1,4)+cb(1,5) . +cb(1,6)+cb(1,7)+cb(1,8) c calculate the derivatives: dCb1(1) = dB1a(1)*asum*exp1/P + B1*dasum*exp1/P . +B1*asum*dexp1/P - B1*asum*exp1*dP(1)/Psq dCb1(2) = dB1a(2)*asum*exp1/P + B1*dasum*exp1/P . +B1*asum*dexp1/P - B1*asum*exp1*dP(2)/Psq dCb1(3) = dB1a(3)*asum*exp1/P + B1*dasum*exp1/P . +B1*asum*dexp1/P - B1*asum*exp1*dP(3)/Psq c dCb2(1) = dbsum*dB1a(1)*exp2 + bsum*dexp2 dCb2(2) = dbsum*dB1a(2)*exp2 + bsum*dexp2 dCb2(3) = dbsum*dB1a(3)*exp2 + bsum*dexp2 c dCb3(1) = dB2(1)*csum + B2*( c31*dB1a(1)*exp1 +c31*B1*dexp1 . +c32 *2.d0*B1*dB1a(1)*exp2 + c32 *B12*dexp2 ) dCb3(2) = dB2(2)*csum + B2*( c31 *dB1a(2)*exp1 +c31 *B1*dexp1 . +c32 *2.d0*B1*dB1a(2)*exp2 + c32 *B12*dexp2 ) dCb3(3) = dB2(3)*csum + B2*( c31 *dB1a(3)*exp1 +c31 *B1*dexp1 . +c32 *2.d0*B1*dB1a(3)*exp2 + c32 *B12*dexp2 ) c dCb4 equations may 27/90: dCb4(1) = dB1a(1)*B3 *dsum1 + B1*dB3(1)*dsum1 + B1*B3 *ddsum1 . + dB1a(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2 dCb4(2) = dB1a(2)*B3 *dsum1 + B1*dB3(2)*dsum1 + B1*B3 *ddsum1 . + dB1a(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2 dCb4(3) = dB1a(3)*B3 *dsum1 + B1*dB3(3)*dsum1 + B1*B3 *ddsum1 . + dB1a(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2 c dCb5 equations may 27/90: (corrected jun01) dCb5(1) = dB1a(1)*fsum*exp7/P + B1*dfsum*exp7/P . + B1*fsum*dexp7/P - B1*fsum*exp7*dP(1)/Psq dCb5(2) = dB1a(2)*fsum*exp7/P + B1*dfsum*exp7/P . + B1*fsum*dexp7/P - B1*Fsum*exp7*dP(2)/Psq dCb5(3) = dB1a(3)*fsum*exp7/P + B1*dfsum*exp7/P . + B1*fsum*dexp7/P - B1*Fsum*exp7*dP(3)/Psq c dCb6(1) = dB1a(1)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7 dCb6(2) = dB1a(2)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7 dCb6(3) = dB1a(3)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7 c dCb7(1) = dB1a(1)*aasum*exp2 + B1*daasum*dP(1)*exp2 . + B1* aasum*dexp2 dCb7(2) = dB1a(2)*aasum*exp2 + B1*daasum*dP(2)*exp2 . + B1* aasum*dexp2 dCb7(3) = dB1a(3)*aasum*exp2 + B1*daasum*dP(3)*exp2 . + B1* aasum*dexp2 dCb8(1) = dB1a(1)*ffsum*exp7 +B1*dffsum*dP(1)*exp7 . +B1* ffsum*dexp7 dCb8(2) = dB1a(2)*ffsum*exp7 +B1*dffsum*dP(2)*exp7 . +B1* ffsum*dexp7 dCb8(3) = dB1a(3)*ffsum*exp7 +B1*dffsum*dP(3)*exp7 . +B1* ffsum*dexp7 dCbnds(1,1) = dCb1(1) +dCb2(1) +dCb3(1) +dCb4(1) +dCb5(1) . +dCb6(1) +dCb7(1) +dCb8(1) dCbnds(1,2) = dCb1(2) +dCb2(2) +dCb3(2) +dCb4(2) +dCb5(2) . +dCb6(2) +dCb7(2) +dCb8(2) dCbnds(1,3) = dCb1(3) +dCb2(3) +dCb3(3) +dCb4(3) +dCb5(3) . +dCb6(3) +dCb7(3) +dCb8(3) c c ----- calculate the Cbndb correction term ---------------- B1 = B1b B12 = B1*B1 B13 = B12*B1 B14 = B13*B1 B15 = B14*B1 asum = d11 + d12*R + d13*Rsq + d14/R + d15/Rsq bsum = d21*B12 + d22*B13 + d23*B14 + d24*B15 csum = d31*B1 *exp1 + d32*B12*exp2 dsum1 = d41*exp1 + d42*exp2 dsum2 = d43*exp1 + d44*exp2 fsum = dx1 + d52*R + d53*Rsq gsum = d61 + d62/R + d63/Rsq aasum = d71 + d72*P + d73*Psq + d74/P + d75/Psq ffsum = d81*P + d82*Psq + d84/Psq dasum = d12 + 2.d0*d13*R -d14/Rsq -2.d0*d15/Rcu dbsum = 2.d0*d21*B1 +3.d0*d22*B12 +4.d0*d23*B13 +5.d0*d24*B14 ddsum1= d41*dexp1 + d42*dexp2 ddsum2= d43*dexp1 + d44*dexp2 dfsum = d52 + 2.d0*d53*R dgsum = -d62/Rsq - 2.d0*d63/Rcu daasum = d72 +2.d0*d73*P -d74/Psq -2.d0*d75/Pcu dffsum = d81 + 2.d0*d82*P -2.d0*d84/Pcu cb(2,1) = B1*asum * exp1 /P cb(2,2) = bsum * exp2 cb(2,3) = B2*csum cb(2,4) = B1*B3*dsum1 + B1*B3b*dsum2 cb(2,5) = B1 * fsum * exp7 /P cb(2,6) = B1 * gsum * exp7 cb(2,7) = B1 * aasum * exp2 cb(2,8) = B1 * ffsum * exp7 Cbnds(2) = cb(2,1)+cb(2,2)+cb(2,3)+cb(2,4)+cb(2,5) . +cb(2,6)+cb(2,7)+cb(2,8) c calculate the derivatives: dCb1(1) = dB1b(1)*asum*exp1/P + B1*dasum*exp1/P . +B1*asum*dexp1/P - B1*asum*exp1*dP(1)/Psq dCb1(2) = dB1b(2)*asum*exp1/P + B1*dasum*exp1/P . +B1*asum*dexp1/P - B1*asum*exp1*dP(2)/Psq dCb1(3) = dB1b(3)*asum*exp1/P + B1*dasum*exp1/P . +B1*asum*dexp1/P - B1*asum*exp1*dP(3)/Psq c dCb2(1) = dbsum*dB1b(1)*exp2 + bsum*dexp2 dCb2(2) = dbsum*dB1b(2)*exp2 + bsum*dexp2 dCb2(3) = dbsum*dB1b(3)*exp2 + bsum*dexp2 c dCb3(1)= dB2(1)*csum + B2*( d31 *dB1b(1)*exp1 +d31 *B1*dexp1 . +d32 *2.d0*B1*dB1b(1)*exp2 + d32 *B12*dexp2 ) dCb3(2)= dB2(2)*csum + B2*( d31 *dB1b(2)*exp1 +d31 *B1*dexp1 . +d32 *2.d0*B1*dB1b(2)*exp2 + d32 *B12*dexp2 ) dCb3(3)= dB2(3)*csum + B2*( d31 *dB1b(3)*exp1 +d31 *B1*dexp1 . +d32 *2.d0*B1*dB1b(3)*exp2 + d32 *B12*dexp2 ) c dCb4 equations may 27/90: dCb4(1) = dB1b(1)*B3 *dsum1 + B1*dB3(1)*dsum1 + B1*B3 *ddsum1 . + dB1b(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2 dCb4(2) = dB1b(2)*B3 *dsum1 + B1*dB3(2)*dsum1 + B1*B3 *ddsum1 . + dB1b(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2 dCb4(3) = dB1b(3)*B3 *dsum1 + B1*dB3(3)*dsum1 + B1*B3 *ddsum1 . + dB1b(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2 c dCb5 equations may 27/90: dCb5(1) = dB1b(1)*fsum*exp7/P + B1*dfsum*exp7/P . + B1*fsum*dexp7/P - B1*fsum*exp7*dP(1)/Psq dCb5(2) = dB1b(2)*fsum*exp7/P + B1*dfsum*exp7/P . + B1*fsum*dexp7/P - B1*fsum*exp7*dP(2)/Psq dCb5(3) = dB1b(3)*fsum*exp7/P + B1*dfsum*exp7/P . + B1*fsum*dexp7/P - B1*fsum*exp7*dP(3)/Psq c dCb6(1) = dB1b(1)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7 dCb6(2) = dB1b(2)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7 dCb6(3) = dB1b(3)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7 c dCb7(1) = dB1b(1)*aasum*exp2 + B1*daasum*dP(1)*exp2 . + B1* aasum*dexp2 dCb7(2) = dB1b(2)*aasum*exp2 + B1*daasum*dP(2)*exp2 . + B1* aasum*dexp2 dCb7(3) = dB1b(3)*aasum*exp2 + B1*daasum*dP(3)*exp2 . + B1* aasum*dexp2 dCb8(1) = dB1b(1)*ffsum*exp7 +B1*dffsum*dP(1)*exp7 . +B1* ffsum*dexp7 dCb8(2) = dB1b(2)*ffsum*exp7 +B1*dffsum*dP(2)*exp7 . +B1* ffsum*dexp7 dCb8(3) = dB1b(3)*ffsum*exp7 +B1*dffsum*dP(3)*exp7 . +B1* ffsum*dexp7 dCbnds(2,1) = dCb1(1) +dCb2(1) +dCb3(1) +dCb4(1) +dCb5(1) . +dCb6(1) +dCb7(1) +dCb8(1) dCbnds(2,2) = dCb1(2) +dCb2(2) +dCb3(2) +dCb4(2) +dCb5(2) . +dCb6(2) +dCb7(2) +dCb8(2) dCbnds(2,3) = dCb1(3) +dCb2(3) +dCb3(3) +dCb4(3) +dCb5(3) . +dCb6(3) +dCb7(3) +dCb8(3) c c calculate the total derivative from the pieces: dCbnda(1) = dT(1) * Cbnds(1) + sumT * dCbnds(1,1) dCbndb(1) = dT(1) * Cbnds(2) + sumT * dCbnds(2,1) dCbnda(2) = dT(2) * Cbnds(1) + sumT * dCbnds(1,2) dCbndb(2) = dT(2) * Cbnds(2) + sumT * dCbnds(2,2) dCbnda(3) = dT(3) * Cbnds(1) + sumT * dCbnds(1,3) dCbndb(3) = dT(3) * Cbnds(2) + sumT * dCbnds(2,3) Cbnda = sumT*Cbnds(1) Cbndb = sumT*Cbnds(2) c some debugging print statements: d if(ipr.gt.1)then d write(7,7650) B2,B3 d write(7,7700) db2,db3 d write(7,*) d write(7,*) 'Vbend values :' d k=1 d write(7,7150) k,dvb1,dvb2,dvb3,dvb4,dvb5,dVbnda d write(7,7155) (vb(1,i),i=1,5) d k=2 d write(7,7150) k,dvb1,dvb2,dvb3,dvb4,dvb5,dVbndb d write(7,7155) (vb(2,i),i=1,5) d write(7,7100) Vbnda,Vbndb d write(7,*) d write(7,*) 'Cbend values:' d write(7,*) 'c81,c81,c82=',c81,c81,c82 d write(7,*) 'c91,c83,c84=',c91,c83,c84 d write(7,*) 'd81,d81,d82=',d81,d81,d82 d write(7,*) 'd91,d83,d84=',d91,d83,d84 d k=1 d write(7,2020) k,(cb(1,j),j=1,9) d write(7,2033) sumT,Cbnds(1) d k=2 d write(7,2020) k,(cb(k,j),j=1,9) d write(7,2033) sumT,Cbnds(2) d k=1 d write(7,7152) k,dcb1,dcb2,dcb3,dcb4,dcb5,dcb6, d . dcb7,dcb8,dcb9,(dCbnds(1,i),i=1,3) d k=2 d write(7,7152) k,dcb1,dcb2,dcb3,dcb4,dcb5,dcb6, d . dcb7,dcb8,dcb9,(dCbnds(k,i),i=1,3) d write(7,7200) dT,dCbnda,dCbndb d write(7,*) 'VbCb_HHH_BKMP exit -------------------------' d write(7,*) d end if return 2000 format(5x,'r1, r2, r3 = ',3(1x,f12.8)) 2010 format(5x,'B1a, B1b, B2, B3 = ',4(1x,g12.6)) 7100 format(' Vbnda = ',g12.6,' Vbndb = ',g12.6) 7650 format(' B2 = ',g12.6,' B3 = ',g12.6) 7700 format(' dB2 = ',3(1x,g12.6),/, . ' dB3 = ',6(1x,g12.6)) 7150 format(' k=',i1,' derivatives: ',/, . ' dVb1 = ',3(1x,g16.10),/, . ' dVb2 = ',3(1x,g16.10),/, . ' dVb3 = ',3(1x,g16.10),/, . ' dVb4 = ',3(1x,g16.10),/, . ' dVb5 = ',3(1x,g16.10),/, . ' dVbnd = ',3(1x,g16.10)) 7155 format(' Vb1 = ',g16.10,' Vb2 = ',g16.10,/, . ' Vb3 = ',g16.10,' Vb4 = ',g16.10,/, . ' Vb5 = ',g16.10) 2020 format(5x,'k=',i1,' cb1,cb2,cb3= ',3(1x,g16.10), . /,8x,' cb4,cb5,cb6= ',3(1x,g16.10), . /,8x,' cb7,cb8,cb9= ',3(1x,g16.10)) 2033 format(5x,'sumT = ',g12.6 ,' cbnds(k) = ',g12.6) 7152 format('Subr.Cbend: k=',i1,' derivatives: ',/, . ' dCb1 = ',3(1x,g18.12),/, . ' dCb2 = ',3(1x,g18.12),/, . ' dCb3 = ',3(1x,g18.12),/, . ' dCb4 = ',3(1x,g18.12),/, . ' dCb5 = ',3(1x,g18.12),/, . ' dCb6 = ',3(1x,g18.12),/, . ' dCb7 = ',3(1x,g18.12),/, . ' dCb8 = ',3(1x,g18.12),/, . ' dCb9 = ',3(1x,g18.12),/, . ' dCbnd = ',3(1x,g18.12)) 7200 format(' dT = ',3(1x,g18.12),/, . ' dCbnda= ',3(1x,g18.12),/, . ' dCbndb= ',3(1x,g18.12)) 6000 format(' from subr.b1ab: r=',3(1x,f9.6)) 6010 format(' b1a = ',f16.10,' b1b = ',f16.10) 6020 format(' db1a: ',3(1x,e12.6),/,' db1b: ',3(1x,e12.6)) 7000 format(' dc(j)/dr(i) derivatives: ',/, . 8x,3(1x,e12.6),/,8x,3(1x,e12.6),/,8x,3(1x,e12.6)) 7010 format(' ri2-rj2-rk2: ',3(1x,e12.6),/, . ' cosines : ',3(1x,e12.6),/, . ' 12c**2 - 3 : ',3(1x,e12.6)) end c