c*********************************************************************** Program betaFIT16 c*********************************************************************** c+++++ COPYRIGHT 2016 by R.J. Le Roy and Asen Pashov ++++++++++++++++++ c++ as per our paper in the 'Remote Sensing' issue of JQSRT (2016) +++++ c+++++++++++++++ version of 25 March 2016 +++++++++++++++++++++++++++++ c* Program to fit NTP read-in potential fx. values {RTP(i),VTP(i)} c (with or without individual weights) to a chosen analytic form. c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cc INCLUDE 'arrsizes.h' c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c** This 'Block' Data Utility routine with the name 'arrsizes.h', that c governs array dimensioning in program betaFIT MUST be installed c under this name in the same (sub)directory containing the FORTRAN c file(s) for this Program when it is being compiled, or incorporated c into the program wherever dimensioning is required (as it is below) c----------------------------------------------------------------------- INTEGER MXDATA, LMAX, MXPARM, NCMMAX PARAMETER (MXDATA= 1001, LMAX= MXDATA, MXPARM=50, NCMMAX= 25) c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c======================================================================= INTEGER i,j,k,INFL,ITER,IROUND,JROUND,LPPOT,ROBUST,prFIT,prNLL, 1 prDIFF,M,NPARM,NPR,NTP,NLIN,NS,NL,NDGF,NxPSE,CYCMAX,IFXP(MXPARM) REAL*8 PV(MXPARM),PU(MXPARM),PS(MXPARM), 1 CM(MXPARM,MXPARM),DYDP(MXDATA,MXPARM),VTP(MXDATA),ypeq(MXDATA), 2 uVTP(MXDATA),betay(MXDATA),Ubetay(MXDATA),YD(MXDATA), 3 yqPSE(MXDATA),xPSE(MXDATA),rPSE(MXDATA),rKL(1:MXDATA,1:MXDATA), 4 DM(NCMMAX),DMP(NCMMAX),DMPP(NCMMAX),ypRegrid(MXDATA), 3 betaINF,UNC,yPOW,DSE,TSTPS,TSTPU,DSEB, TCM,UM,diff, 4 Rep,Req,AREF,AREFp,AREFq,RTPp,RTPq, ULR,FCT,RAT,UMAX, 5 XX,YY,YH,yp,yq,ypRE,ReDE, ReIN,DeIN,VMINin,ULRe,dULRedRe, 6 dULRdR,RE3,RE6,RE8,T0,T1,C6adj,C9adj,RTP3,RTP6,RTP8,RH,RR,RB,RBB, 7 VV,VB,VBB,SCALC,yMIN,DRMSD ,BETA0,BETAN, RPR1,dRPR, SL,SLB,dSL, 8 ycalc,dd,VLIM,yqREFre REAL*8 DEIGM1(1,1),DEIGM3(1,1),DEIGM5(1,1),DEIGR(1,1),DEIGRe(1,1) 1 ,DEIGDe(1,1) CHARACTER*4 NNAME,NAME(4) DATA NAME/' EMO',' MLR','DELR','GPEF'/,CYCMAX/25/ c----------------------------------------------------------------------- REAL*8 Re,De,VMIN,RREF,AA,BB,as,bs,rhoAB,BETA, CmVAL(NCMMAX), 1 CmEFF(NCMMAX),dULRdCm(NCMMAX),RTP(MXDATA),SAS(MXDATA,MXPARM) INTEGER PSEL,IFXRe,IFXDe,IFXVMIN,sVSR2,IDSTT,NCMM,MCMM,q,p, 1 NPbeta,APSE,MMLR(NCMMAX),IFXCm(NCMMAX),updateCmADJ COMMON /DATABLK/Re,De,VMIN,RREF,AA,BB,as,bs,rhoAB,CmVAL,CmEFF, 1 dULRdCm,RTP,SAS,BETA,PSEL,IFXRe,IFXDe,IFXVMIN,sVSR2,IDSTT,NCMM, 2 MCMM,q,p,NPbeta,APSE,MMLR,IFXCm,updateCmADJ c----------------------------------------------------------------------- ROBUST= 0 c----------------------------------------------------------------------- c** PSEL specifies the type of potential being fitted to: c PSEL= 1 for EMO; PSEL=2 for an MLR (or MLJ); PSEL=3 for DELR ; c PSEL= 4 for GPEF c* NPT is the number of read-in potential points being fitted to. c* UNC is the energy uncertainty associated with the potential points c (plausibly ca. 0.1 cm-1 for RKR). To weight each point separately c set UNC < 0.0 and read in a separate uncertainty for each point c IROUND specifies the level of rounding inside NLLSSRR if: c > 0 : requires that Sequential Rounding & Refitting be c performed, with each parameter being rounded at the c IROUND'th sig. digit of its local uncertainty. c <=0 : simply stops after full convergence (without rounding). c LPPOT specifies whether (LPPOT > 0) or not (LPPOT.le.0) the code should c read instructions to generate an array of potential function c values on some specified mesh and range c If LPPOT > 0, read a separate mesh and range for eash {N,q,p,Rref} case c* prFIT = 1 causes printout of results of preliminary linearized fit c and (as appropriate) other non-final fits. Normally set prFIT=0 c prFIT > 1 specifies the level of printing inside fit subroutine NLLSSRR c >= 1 print converged parameters, PU & PS's c >= 2 also print parameter change each rounding step c >= 3 also indicate nature of convergence c >= 4 also print convergence tests on each cycle c >= 5 also parameters changes & uncertainties, each cycle c prDIFF (integer) > 0 causes fit residuals to be printed after each c complete fit: no print if prDIFF.LE.0 c** Re & De are assumed potential well minimum position & well depth. c ... note that De is a dummy variable for a GPEF function c** VMIN is the minimum of the potential function defined by the read-in c points [typically =0 for RKR potential, but non-zero for ab initio] c** To fix Re, De or VMIN unchanged at read-in values, set (integer) c IFXRe, IFXDe, and/or IFXVMIN > 0 ; to fit them, set value =0 c [Normally set IFXVMIN= 0 !!] c======================================================================= READ(5,*) PSEL, NTP, UNC, IROUND, LPPOT, prFIT, prDIFF READ(5,*) Re, De, VMIN READ(5,*) IFXRe, IFXDe, IFXVMIN c======================================================================= JROUND= 0 prNLL= 0 IF(prFIT.ge.2) prNLL= 5 NNAME= NAME(PSEL) ReIN= Re DeIN= De VMINin= VMIN IF(PSEL.LE.3) WRITE(6,600) NNAME, VMIN, Re, De IF(PSEL.GT.3) WRITE(6,600) NNAME, VMIN, Re c** For an MLR (PSEL=2) or DELR (PSEL=3) potential, read the number of c long-range terms NCMM used to define long-range potential tail: c V(r)= De - \sum{CmVAL/r^MMLR} c** If rhoAB .LE. 0.0 have NO damping functions: all Dm(R)= 1.0 c If rhoAB > 0.0 it is the molecule-dependent radial scaling factor c of Douketis et al. [JCP 76, 3057 (1982)] c rhoAB = 2*rhoA*rhoB/(rhoA+rhoB) where rhoA is the ionization c potential ratio (I_p^A/I_p^H)^{2/3} for atom A vs. atomic H c For rhoAB > 0.0, sVSR2 specifies damping s.th. Dm(r)/r^m --> r^{sVSR2/2} c IDSTT > 0 use generalized Douketis et al. damping fx. c IDSTT > 1 use DS damping for s=-1/2 AND impose very c short-range 1/r behaviour with C1= IDSTT = Z1*Z2 c IDSTT .LE. 0 use generalized Tang-Toennies damping fx c** APSE is an integer: > 0 to use A.Pashov natural Spline in MLR Exponent c : .le. 0 for normal constrained polynomial exponent c* For an MLR potential with APSE > 0, yMIN is a negative{!} number c (-1 .le. yMIN < 0) which specifies the lower bound on the yp values c selected to define the MLR exponent spline. c** For each long-range term read power MMLR(i) & coefficient CmVAL(i) c** For special Aubert-Frecon 2x2 cases, NCMM= 7, MMLR= {x,3,3,6,6,8,8}, c with x= 0 for the A state, x= -1 for the b state, and CmVAL= {Aso, c C3Sig, C3Pi, C6Sig, C6Pi, C8Sig, C8Pi}, c* while for the 3x3 diagonalization cases, NCMM=10, MMLR= {x,3,3,3,6,6,6, c 8,8,8} with x= -2 for the (lowest eigenvalue) c(1\,^3\Sigma_g^+ state, c x= -3 for the (middle root) B^1\Pi_u state, and x=-4 for the c highest-root state, while CmVal= {Aso, C3Sig, C3Pi1, C3Pi3, C6Sig, C6Pi1, c C6Pi3, C8Sig, C8Pi1, C8Pi3} c======================================================================= IF((PSEL.EQ.2).OR.(PSEL.EQ.3)) THEN READ(5,*) NCMM, rhoAB, sVSR2, IDSTT, APSE, yMIN DO m= 1, NCMM READ(5,*) MMLR(m),CmVAL(m) CmEFF(m)= CmVAL(m) ENDDO ENDIF MCMM= NCMM IF(PSEL.NE.2) APSE= 0 c======================================================================= c** For a GPEF potential, read coefficients to define expansion vble: c y = (r^q - Re^q)/(as*r^q + bs*Re^q) where q, as & bs all fixed c======================================================================= IF(PSEL.EQ.4) THEN READ(5,*) as, bs WRITE(6,604) as,bs ENDIF c======================================================================= c** Read the turning points to be fitted to c======================================================================= IF(UNC.GT.0.d0) READ(5,*) (RTP(i), VTP(i),i= 1,NTP) IF(UNC.LE.0.d0) READ(5,*) (RTP(i), VTP(i),uVTP(i),i= 1,NTP) c======================================================================= IF(UNC.GT.0.d0) THEN WRITE(6,606) NTP,UNC,(RTP(i),VTP(i),i= 1,NTP) DO i= 1,NTP uVTP(i)= UNC ENDDO ELSE WRITE(6,609) NTP,VMIN,(RTP(i),VTP(i),uVTP(i),i= 1,NTP) ENDIF WRITE(6,608) IF((NCMM.GE.4).AND.((MMLR(1).EQ.0).OR.(MMLR(1).EQ.-1))) THEN c** Shift of asymptote for special alkali dimer 2x2 cases DO I= 1,NTP !! ???? huh?? VTP(I)= VTP(I) - 0.5d0*CmVAL(1) !! ???? huh?? ENDDO !! ???? huh?? ENDIF c======================================================================= c 690 FORMAT(' uLR(r) is a damped TT inverse-power sum with',6x, 1' rhoAB=', F7.4,' sVSR2=',i3/(49x,'C',i2,'=',1PD15.8:)) 695 FORMAT(' uLR(r) inverse-power terms incorporate ',a2,' damping wit 1h rhoAB=',f10.7/8x ,'defined to give very short-range damped uLR- 2term behaviour r^{',i2,'/2}') 696 FORMAT(' uLR(r) inverse-power terms incorporate ',a2,' damping wit 1h rhoAB=',f10.7/8x ,'defined to give very short-range damped uLR- 2term behaviour r^{',i2,'/2}'/(49x,'C',i2,'=',1PD15.8:)) 698 FORMAT(' uLR(r) is a simple inverse-power sum with',6x,'C',I2, 1 '=',1PD15.8:/(49x,'C',i2,'=',D15.8:)) 600 FORMAT(' Fit an ',A4,' potential function to the input points'/ 1 1x,57('=')/5x,'with initial VMIN=',f13.4,' Re=',f11.8: 2 ,' De=',f11.4) 614 FORMAT(' Fit an ',A4,'(q=',I2,' p=',I2,') potential function to t 1he input points'/1x,57('=')/5x,'with initial VMIN=',f13.4, 2 ' Re=',f11.8:,' De=',f11.4) 615 FORMAT(' Fit an ',A4,'(q=',I2,') potential function to the input p 1oints'/1x,57('=')/5x,'with initial VMIN=',f13.4,' Re=',f11.8: 2 ,' De=',f11.4) 601 FORMAT(6x,'in which y_{p/q}(r)= [r^{p/q} -',f7.4,'^{p/q}]/[r^{p 1/q} +',f7.4,'^{p/q}]' ) 602 FORMAT(5x,'Use Lyon 2x2 ',A7,' uLR(r) with Aso=',F11.6/ 1 47x,'C_3(^1Sig)=',1PD15.7:/47x,'C_3(^3Pi) =',D15.7:/ 1 47x,'C_6(^1Sig)=',1PD15.7:/47x,'C_6(^3Pi) =',D15.7:/ 1 47x,'C_8(^1Sig)=',1PD15.7:/47x,'C_8(^3Pi) =',D15.7) 6022 FORMAT(' Use Lyon 3x3 ',A7,' uLR(r) with Aso=',F11.6 / 1 47x,'C_3(^3Sig)=',D15.7:/47x,'C_3(^1Pi) =',D15.7:/ 2 47x,'C_3(^3Pi) =',D15.7:/ 3 47x,'C_6(^3Sig)=',D15.7:/47x,'C_6(^1Pi) =',D15.7:/ 4 47x,'C_6(^3Pi) =',D15.7:/ 5 47x,'C_8(^3Sig)=',D15.7:/47x,'C_8(^1Pi) =',D15.7:/ 6 47x,'C_8(^3Pi) =',D15.7) 603 FORMAT(6x, 'in which y_{p/q}(r)= [r^{p/q} - Re^{p/q}]/[r^{p/q} + 1 Re^{p/q}]' ) 652 FORMAT( ' & define beta(y_q^{ref}(r)) as a natural spline through 1 points at the',i4,' yq^{ref} values:'/(2x,7F11.7)) 604 FORMAT(' GPEF expansion variable is: y= (R^q - Re^q)//(',F9.5, 1 '*R^q',1x,SP,F9.5,'*Re^q)' ) 605 FORMAT(' GPEF expansion variable is: y= (R^',I1,' - Re^',i1, 1 ')/(',F9.5,'*R^',I1,1x,SP,F9.5,'*Re^',SS,i1,')' ) 606 FORMAT(/' Fit to',I5,' input turning points assuming common uncert 1ainty u(VTP)=',1PD9.2/1x,39('--')/ 2 4(' RTP',7x,'VTP ')/1x,39('--')/(4(0pF9.5,f11.3))) 607 FORMAT(' which yields initial values of AA=',1PD14.7, 1 ' BB=',D14.7) 608 FORMAT(1x,39('--')/) 609 FORMAT(/' Fit to',I5,' input turning points with initial energy mi 1nimum VMIN=',f11.4/1x,32('--')/2(5x,'RTP',8x,'VTP',6x,'unc',4x)/ 2 1x,32('--')/(2(0PF10.5,F12.4,1PD10.2))) 610 FORMAT(' Shift read-in VMIN value below lowest input potential val 1ue',f12.4) 611 FORMAT(3x,'Start with this long-range tail and beta(0)=',f10.6) 612 FORMAT(' *** Linearization Fails ** for datum',I4,' r =',f8.2/ 1 ' De too small or VMIN too high') c======================================================================= c** Now ... loop over different {q,p,NS,NL} combinations till end of data c** p and q are powers used to define radial variable in the exponent c beta(r)= yp*betaINF + [1 - yp]*sum{beta_i*yq^i} where c ya=(R^a - AREF^a)/(R^a + AREF^a) for a= p or q c** NL= is the order of the\beta(y) exponent expansion for EMO, DELR, c or MLR potentials with APSE.leq.0. No. such param is NPbeta= NL+ 1 c NS is a dummy variable for EMO, DELR & 'ordinary' MLR with APSE.le.0 c** For an MLR potential with APSE > 0, NPbeta=(NS+NL+1) is the number c of yp values to be used to define the exponent spline used for beta(y): c NL specifies the number of yp values from yp= 0 to yp= 1.0 c NS specifies the number of (equally spaced) yp values for c yMIN .leq. yp = yMIN to 0.0 c** For a GPEF potential, consider powers ranging from NS(.ge.0) to NL c* RREF defines the reference distance in the expansion variable c - for RREF.le.0 , define parameter RREF = Re c - for RREF.gt.0 , fix parameter RREF at its read-in value c----------------------------------------------------------------------- 20 READ(5,*,END= 999) q, p, NS, NL, RREF c----------------------------------------------------------------------- IF(q.LE.0) GO TO 999 IF(PSEL.NE.2) p = q NPR= 0 c======================================================================= c** If desired, print the resulting potential for this case to channel-8 c at for NPR distances, starting at RPR1 with a mesh of DPR c* IF(NPR.leq.0) omit printout for this case c----------------------------------------------------------------------- IF(LPPOT.GT.0) READ(5,*) NPR, RPR1, dRPR c----------------------------------------------------------------------- Re= ReIN De= DeIN VMIN= VMINin IF(PSEL.LE.3) NPbeta= NL + 1 IF(PSEL.EQ.2) WRITE(6,614) NNAME, q, p, VMIN, Re, De IF((PSEL.EQ.1).OR.(PSEL.EQ.3)) WRITE(6,615) NNAME, q, VMIN, Re, De IF(PSEL.EQ.4) THEN WRITE(6,600) NNAME, VMIN, Re WRITE(6,605) q,q,as,q,bs,q ENDIF IF(PSEL.LE.3) THEN IF(RREF.GT.0.d0) THEN AREF= RREF WRITE(6,632) q,q,RREF,q,q,RREF,q ELSE AREF= RE IF(RREF.LE.0.d0) WRITE(6,633) q,q,q,q,q ENDIF ENDIF AREFp= AREF**p AREFq= AREF**q IF(PSEL.EQ.2) THEN IF(APSE.GT.0) WRITE(6,650) NS, NL IF(APSE.LE.0) WRITE(6,634) p, p, q IF(RREF.gt.0.d0) THEN AREF= RREF WRITE(6,601) AREF,AREF ELSE AREF= Re WRITE(6,603) ENDIF IF(APSE.GT.0) THEN c** Set up yp abscissa values for APSpline Exponent case c** first ... FIND y_q^{ref}(r_e) yqREFre= (RE**q - AREFq)/(RE**q + AREFq) YH= (yqREFre -yMIN)/NS yqPSE(1)= yMIN c** now .. select NS-1 points between yqREFre and YMIN DO I= 2,NS yqPSE(I)= yqPSE(I-1) + YH ENDDO NPbeta= NS+ 1 yqPSE(NS+1)= yqREFre - 0.05d0 c** now .. select NL-1 points from yqREFre to YMAX= 1.0 YH= (1.d0 - yqREFre)/NL DO I= 1,NL-1 yqPSE(NPbeta+I)= yqPSE(NPbeta+I-1) + YH ENDDO NPbeta= NPbeta+ NL yqPSE(NPbeta)= 1.d0 WRITE(6,652) NPbeta,(yqPSE(i), i= 1,NPbeta) DO i= 1,NPbeta c*** Round yp values for 'portability' & to avoid precisely Re J= 100* yqPSE(i) yqPSE(i)= FLOAT(J)*1.d-2 ENDDO WRITE(6,652) NPbeta,(yqPSE(i), i= 1,NPbeta) ENDIF ENDIF IF((NCMM.GE.4).AND.(MMLR(1).LE.0)) THEN c** uLR printout for Lyon 2x2 or 3x3 treatment of 2S + 2p alkali dimers ... IF(MMLR(1).EQ.0) WRITE(6,602) 'A-state',CmVAL(1), 1 CmVAL(2),(CmVAL(m),m=3,NCMM) c... For Lyon treatment of b-state alkali dimers ... IF(MMLR(1).EQ.-1) WRITE(6,602) 'b-state',CmVAL(1), 1 CmVAL(2),(CmVAL(m),m=3,NCMM) IF(MMLR(1).EQ.-2) WRITE(6,6022) 'c-state',CmVAL(1) ,CmVAL(2), 1 (CmVAL(m),m=3,NCMM) IF(MMLR(1).EQ.-3) WRITE(6,6022) 'B-state',CmVAL(1) ,CmVAL(2), 1 (CmVAL(m),m=3,NCMM) IF(MMLR(1).EQ.-4) WRITE(6,6022) '1 ^3Pi',CmVAL(1) ,CmVAL(2), 1 (CmVAL(m),m=3,NCMM) IF(IDSTT.GT.0) WRITE(6,695) 'DS',rhoAB,sVSR2 IF(IDSTT.LE.0) WRITE(6,695) 'TT',rhoAB,sVSR2 ENDIF c... uLR printout for 'conventional' (damped or non-damped) inverse-power sum IF(((PSEL.EQ.3).OR.((PSEL.EQ.2))).AND.(MMLR(1).GT.0)) THEN IF(rhoAB.LE.0) THEN WRITE(6,698) (MMLR(m),CmVAL(m),m= 1,NCMM) ELSE IF(IDSTT.GT.0) WRITE(6,696) 'DS',rhoAB,sVSR2, 1 (MMLR(m),CmVAL(m),m= 1,NCMM) IF(IDSTT.LE.0) WRITE(6,696) 'TT',rhoAB,sVSR2, 1 (MMLR(m),CmVAL(m),m= 1,NCMM) ENDIF ENDIF c======================================================================= IF(PSEL.EQ.2) CALL quadCORR(NCMM,MCMM,NCMMAX,MMLR,De,CmVAL,CmEFF) c======================================================================= DSE= VMIN c** Scan input potential array to ensure VMIN .le. {lowest input point} DO i= 1,NTP DSE= DMIN1(DSE,VTP(i)) ENDDO IF(DSE.LT.VMIN) THEN VMIN= DSE - 0.1d0 WRITE(6,610) VMIN ENDIF c======================================================================= c** Preliminary linearized fit for an EMOp potential ... c----------------------------------------------------------------------- IF(PSEL.EQ.1) THEN c ... first define ordinate array DO i= 1,NTP RTPp= RTP(i)**p yp= (RTPp - AREFp)/(RTPp + AREFp) IF(RTP(i).GT.Re) THEN T0= (1.d0 - DSQRT((VTP(i)-VMIN)/De)) IF(T0.LT.0.d0) THEN WRITE(6,612) i,RTP(i) STOP ENDIF betay(i)= - DLOG(T0) IF(VTP(i).GT.uVTP(i)) THEN Ubetay(i)= 0.5d0*uVTP(i)/(DSQRT((VTP(i)-VMIN)*De) 1 - (VTP(i)-VMIN)) ELSE Ubetay(i)= DSQRT(uVTP(i)/De) ENDIF ELSE betay(i)= - DLOG(1.d0 + DSQRT((VTP(i)-VMIN)/De)) IF(VTP(i).GT.uVTP(i)) THEN Ubetay(i)= 0.5d0*uVTP(i)/(DSQRT((VTP(i)-VMIN)*De) 1 + (VTP(i)-VMIN)) ELSE Ubetay(i)= DSQRT(uVTP(i)/De) ENDIF ENDIF c ... next create partial derivative array for linearized fit ... yPOW= (RTP(i)- Re) DO j= 1, NPbeta DYDP(i,j)= yPOW yPOW= yPOW*yp ENDDO c%% elective printout for testing ccc if(i.eq.1) write(8,700) ccc write(8,702) rtp(i),yp,ypRE,vtp(i),betay(i), ccc 1 (dydp(i,j),j=1,NPbeta) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ENDDO CALL LLSQF(NTP,NPbeta,MXDATA,MXPARM,betay,Ubetay,DYDP,YD,PV, 1 PU,PS,CM,DSE) IF(prFIT.GT.0) WRITE(6,620) NNAME,q,NL,DSE, 1 ('beta',j-1,PV(j),PU(j),PS(j),j= 1,NPbeta) ENDIF c======================================================================= c*** Preliminary linearized fit for an MLR potential ... c----------------------------------------------------------------------- IF(PSEL.EQ.2) THEN IF((MMLR(1).GT.0).AND.(p.LE.(MMLR(NCMM)-MMLR(1)))) 1 WRITE(6,616) p,NCMM,MMLR(NCMM),1,MMLR(1) IF((MMLR(1).LE.0).AND.(p.LE.(MMLR(NCMM)-MMLR(2)))) 1 WRITE(6,616) p,NCMM,MMLR(NCMM),2,MMLR(2) c** First define array of exponent values with uncertainties defined by c the assumption that all potential values have equal uncertainties UNC c... Begin by determining uLR(Re) and betaINF IF((NCMM.GE.4).AND.(MMLR(1).LE.0)) THEN c** For special Aubert-Frecon 2x2 or 3x3 alkali dimer cases ... VLIM= VMIN+ De CALL AFdiag(Re,VLIM,NCMM,NCMMAX,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULRe,dULRdCm,dULRedRe) ELSE c... for normal MLR/MLJ case ... with or without damping --------------- CALL dampF(Re,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,DM,DMP, 1 DMPP) ULRe= 0.d0 DO m= 1,NCMM ULRe= ULRe + DM(m)*CmEFF(m)/Re**MMLR(m) ENDDO ENDIF c*** Now, generate betaINF from the adjusted constants ===== betaINF= DLOG(2.d0*De/ULRe) WRITE(6,619) betaINF c.... Have Completed determination of uLR(Re) and betaINF ============ c.... Now loop to prepare linearized MLR arguments & derivativesa------ Rep= RE**p k= 0 DO i= 1, NTP RTPp= RTP(i)**p RTPq= RTP(i)**q yp= (RTPp - AREFp)/(RTPp + AREFp) yq= (RTPq - AREFq)/(RTPq + AREFq) ypRE= (RTPp - Rep)/(RTPp + Rep) ypeq(i)= ypRE cc IF(DABS(yq).le.0.1d0) GO TO 66 !?? omit point? k=k+ 1 xPSE(k)= yq !!placed after the k counter IF((NCMM.GE.4).AND.(MMLR(1).LE.0)) THEN c ... for Aubert-Frecon Li2(A) 2x2 or 3x3 cases CALL AFdiag(RTP(i),VLIM,NCMM,NCMMax,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULR,dULRdCm,dULRdR) ELSE c... for normal inverse-power-sum MLR/MLJ case ... with or without damping CALL dampF(RTP(i),rhoAB,NCMM,NCMMAX,MMLR, 1 sVSR2,IDSTT,DM,DMP,DMPP) ULR= 0.d0 DO m= 1,NCMM ULR= ULR + DM(m)*CmEFF(m)/RTP(i)**MMLR(m) ENDDO ENDIF RAT= DSQRT((VTP(i)-VMIN)/De) IF(RTP(i).GT.Re) THEN T0= (1.d0 - RAT) IF(T0.LT.0.d0) THEN WRITE(6,612) i,RTP(i) STOP ENDIF betay(k)= - DLOG(T0*ULRe/ULR) !! = beta*yp{eq} Ubetay(k)= 0.5d0*uVTP(i)/(RAT*(1d0 - RAT)) ELSE betay(k)= - DLOG((1.d0 + RAT)*ULRe/ULR) Ubetay(k)= 0.5d0*uVTP(i)/(RAT*(1.d0+ RAT)) ENDIF IF(APSE.LE.0) THEN c** Subtract the \beta_\infty term to yield polynomial for fitting c... For Huang MLR polynomial exponent function betay(k)= betay(k)- betaINF*yp*ypRE yPOW= ypRE*(1.d0- yp) c... then create partial derivative array for linearized fit ... DO j= 1, NPbeta DYDP(i,j)= yPOW yPOW= yPOW*yq ENDDO ENDIF c%%%% printout for testing polynomial exponent %%%%%%%%%%%%%%%%%%%%%%%% c if(i.eq.1) write(7,700) c 700 format(' RTP yp yp^{eq} VTP(i) uLR', c 1 10x,' beta(i)') c write(7,702) rtp(i),yp,ypRE,vtp(i),ULR,betay(i), c 1 (dydp(i,j),j=1,NPbeta) c 702 format(f6.3,2f9.5,1P,3d13.5:/(14x,5d13.5)) c%%%% END of printout for testing polynomial exponent %%%%%%%%%%%%%%%%% 66 CONTINUE ENDDO IF(APSE.LE.0) THEN CALL LLSQF(NTP,NPbeta,MXDATA,MXPARM,betay,Ubetay,DYDP,YD, 1 PV,PU,PS,CM,DSE) IF(prFIT.GT.0) THEN IF(RREF.GE.0.d0) WRITE(6,621) NNAME,q,p,RREF,NL, 1 DSE,('beta',j-1,PV(j),PU(j),PS(j),j= 1,NPbeta) IF(RREF.LT.0.d0) WRITE(6,623) NNAME,q,p,NL, 1 DSE,('beta',j-1,PV(j),PU(j),PS(j),j= 1,NPbeta) ENDIF ENDIF IF(APSE.GT.0) THEN c** For Pashov spline-exponent, use spline through linearized input c exponent values to define initial trial spline values of \beta(i) c at the specified xPSE(i) 0:qvalues ! c... First, add xPSE(NTP+1)= 1.0 where betay= \betaINF array point NxPSE= k + 1 xPSE(NxPSE)= 1.d0 betay(NxPSE)= betaINF ubetay(NxPSE)= 1.d0 rtp(NxPSE)= 9.d99 ypeq(NxPSE)= 1.d0 c** Elective printout for linearized SE-MLR stuff ********************* cc write(8,699) (rtp(i),xPSE(i),betay(i),betay(i)/ypeq(i), cc 1 Ubetay(i),ypeq(i), I= 1,NxPSE) cc699 FORMAT(' R=', 0p,f10.6,' xPSE=',F12.9,' betay=', f12.9, cc 1 ' beta=',f15.9,' ubetay=',1pd10.3,' ypeq=',0p,f12.9) c****** end of elective printout *************************************** CALL Lkoef(NxPSE,xPSE,rKL) DO I= 1, NPbeta- 1 !! skip endpoint where rPSE= \infty XX= yqPSE(I) rPSE(I)= AREF*((1.d0 + XX)/ 1 (1.d0 - XX))**(1.d0/DFLOAT(q)) c... Now, use a spline through the exponent values defined by the input c points to generate values of that exponent at the desired c spline-definition points ypREgrid(I)= (rPSE(I)**p - Re**p)/(rPSE(I)**p+ Re**p) PV(I)= 0.d0 DO m= 1, NxPSE PV(I)= PV(I) + 1 Scalc(XX,m,NxPSE,xPSE,rKL,MXDATA) * betay(m) ENDDO ENDDO rPSE(NPbeta)= 9.d99 ypREgrid(NPbeta)= 1.d0 PV(NPbeta)= betaINF !! now include known endpoint!! IF(prFIT.GT.0) WRITE(6,653) NNAME,q,p, 1 (yqPSE(I),rPSE(I),PV(I),ypREgrid(i),PV(I)/ypREgrid(I), 1 I= 1,NPbeta) c*********** Finally, convert beta*ypRE to beta DO I= 1, NPbeta PV(I)= PV(I)/ypREgrid(i) ENDDO c c... Finally, create the fixed array of S(m,x) to define the exponent c and its partial derivatives in subsequent fits, CALL Lkoef(NPbeta,yqPSE,rKL) c** Now Compare predictions of our preliminary spline function defined by c these NPbeta points with the input data ... dd= 0.d0 DO I= 1,NxPSE ycalc= 0.d0 DO m= 1,NPbeta SAS(I,m)= Scalc(xPSE(I),m,NPbeta,yqPSE,rKL,MXDATA) ycalc= ycalc+ SAS(I,m)*PV(m) ENDDO c?? 12/01/15 hey ... maybe need to divide ycalc by xPSE ycalc= ycalc*xPSE(I) !! try this RJL ! c?? It appears to work ... NOW Y_{calc} \approx Y_{input} c?? .. but is this beta or beta*yp or ...?? c?? Now need to fix up ubetay(i) !!! diff= (ycalc - betay(I))/ubetay(I) dd= dd+ diff**2 cc WRITE(6,655) xPSE(I),ycalc,betay(I)/xPSE(I),diff ENDDO dd= DSQRT(dd/NxPSE) WRITE(6,656) NS, NL, RREF, dd ENDIF ENDIF cc655 FORMAT(' At X=',F8.4,' Y_{calc}=',F11.7,' Y_{input}=', cc 1 F11.7,1x,F13.8) 656 FORMAT(' SE-MLR Linearization: NS=',I2,', NL=',I3, 1 ', R_{ref}=', F7.3,' yields dd=',F9.5/) c======================================================================= c** Preliminary linearized fit for a DELR potential ... c----------------------------------------------------------------------- IF(PSEL.EQ.3) THEN BETA0= 1.d0 WRITE(6,611) BETA0 IF((NCMM.GE.4).AND.(MMLR(1).LE.0)) THEN c** For special Aubert-Frecon 2x2 or 3x3 alkali dimer cases ... VLIM= VMIN+ De CALL AFdiag(Re,VLIM,NCMM,NCMMAX,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULRe,dULRdCm,dULRedRe) ELSE CALL dampF(Re,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,DM,DMP, 1 DMPP) ULRe= 0.d0 dULRedRe= 0.d0 DO m= 1,NCMM ULRe= ULRe+ AA*DM(m) dULRedRe= dULRedRe+ AA*(DMP(m) - MMLR(m)*DM(m)/Re) ENDDO ENDIF c ... NOTE need iteration to determine self-consistent beta(0) value c First generate A & B from input Re, De and trial beta(0) 40 CONTINUE AA= De - ULRe - dULRedRe/beta0 BB= 2.d0*(De - ULRe) - dULRedRe/beta0 WRITE(6,607) AA,BB RAT= 0.5d0*BB/AA UMAX= DSQRT(RAT**2 + (uVTP(i) + ULRe - DE)/AA) ReDE= Re- dlog(RAT)/beta0 DO i= 1,NTP ULR= 0.d0 RTPp= RTP(i)**p yp= (RTPp - AREFp)/(RTPp + AREFp) IF((NCMM.GE.4).AND.(MMLR(1).LE.0)) THEN c** For special Aubert-Frecon 2x2 or 3x3 alkali dimer cases ... VLIM= VMIN+ De CALL AFdiag(RTP(i),VLIM,NCMM,NCMMAX,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULR,dULRdCm,dULRdR) ELSE CALL dampF(RTP(i),rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT, 1 DM,DMP,DMPP) DO m= 1,NCMM ULR= ULR + DM(m)*CmEFF(m)/RTP(i)**MMLR(m) ENDDO ENDIF FCT= (VTP(i) - VMIN + ULR - De)/AA + RAT**2 IF(FCT.LT.0.d0) THEN c** If estimate of Re/DE off a bit and FCT < 0 , ignore & deweight point betay(i)= 0.d0 Ubetay(i)= 9.d99 GOTO 44 ENDIF FCT= DSQRT(FCT) IF(RTP(i).GT.ReDE) THEN IF(RAT.GT.FCT) THEN betay(i)= - DLOG(RAT - FCT) IF((VTP(i)-VMIN).GT.uVTP(i)) THEN Ubetay(i)= 0.5d0*uVTP(i)/(AA*(RAT- FCT)*FCT) ELSE Ubetay(i)= uVTP(i)/(AA*UMAX*(RAT- UMAX)) ENDIF ELSE c ... deweight away points for which \ln argument would be negative betay(i)= 0.d0 Ubetay(i)= 9.d9 ENDIF ELSE betay(i)= - DLOG(RAT + FCT) IF((VTP(i)-VMIN).GT.uVTP(i)) THEN Ubetay(i)= 0.5d0*uVTP(i)/(AA*(RAT+ FCT)*FCT) ELSE Ubetay(i)= uVTP(i)/(AA*UMAX*(RAT+ UMAX)) ENDIF ENDIF c ... now create partial derivative array for linearized fit ... 44 yPOW= (RTP(i)- Re) DO j= 1, NPbeta DYDP(i,j)= yPOW yPOW= yPOW*yp ENDDO ENDDO CALL LLSQF(NTP,NPbeta,MXDATA,MXPARM,betay,Ubetay,DYDP,YD,PV, 1 PU,PS,CM,DSE) IF(DABS(PV(1)-BETA0).GT.PS(1)) THEN WRITE(6,644) beta0,PV(1),PV(1)-beta0,DSE beta0= PV(1) ITER= ITER+ 1 IF(ITER.LE.20) GOTO 40 WRITE(6,646) ITER ELSE WRITE(6,648) PV(1),PV(1)-beta0 ENDIF IF(prFIT.GT.0) WRITE(6,620) NNAME,q,NL,DSE, 1 ('beta',j-1,PV(j),PU(j),PS(j),j= 1,NPbeta) ENDIF IF(PSEL.LE.3) THEN c====================================================================== c** Now ... do direct non-linear fits to potential values ... first c with Re and/or VMIN and De fixed, and then freeing them up too ... c======================================================================= c* FIRST optimize BETA(j)'s (and VMIN) with Re and De held fixed! DO j= 1,NPbeta IFXP(j)= 0 ENDDO c**?? Is this how I implement fixed limiting C1 for MLR? IF(IDSTT.GT.1) IFXP(NPbeta)= IDSTT c** Fix parameter value for yp= 1 in SE-MLR IF(APSE.GT.0) IFXP(NPbeta)= 1 c** For initial run, fix the VMIN, De and Re values NPARM= NPbeta+ 3 IFXP(NPARM-2)= 1 IFXP(NPARM-1)= 1 IFXP(NPARM)= 1 NDGF= NTP- NPARM IF(IFXVMIN.LE.0) THEN IFXP(NPARM)= 0 NDGF= NDGF-1 ENDIF PV(NPARM-2)= Re PV(NPARM-1)= De PV(NPARM)= VMIN updateCmADJ= 1 c ..... On first call to NLLSSRR, free only VMIN and the \beta_i CALL NLLSSRR(NTP,NPARM,MXPARM,CYCMAX,JROUND,ROBUST,prNLL,IFXP, 1 VTP,uVTP,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) DRMSD= DSE*DSQRT(DFLOAT(NDGF)/DFLOAT(NTP)) VMIN= PV(NPARM) IF(prFIT.GT.0) THEN IF(APSE.LE.0) THEN IF(RREF.GT.0.d0) THEN WRITE(6,622) NNAME,q,p,AREF,NL,DRMSD,0,PV(1), 1 PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) ENDIF IF(RREF.LE.0.d0) THEN WRITE(6,624) NNAME,q,p,NL,DRMSD,0,PV(1),PU(1), 1 PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) ENDIF ELSE IF(RREF.GT.0.d0) WRITE(6,626) NNAME,q,AREF,NS, 1 NL,DRMSD,(j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) IF(RREF.LE.0.d0) WRITE(6,628) NNAME,q,NS,NL, 1 DRMSD,(j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) ENDIF IF(IFXVMIN.LE.0) THEN WRITE(6,660) PV(NPARM),PU(NPARM),PS(NPARM) ENDIF ENDIF c ... then, if appropriate, repeat NLLSSRR call with Re free too ... IF(IFXRe.LE.0) THEN NDGF= NDGF- 1 IFXP(NPARM-2)= 0 CALL NLLSSRR(NTP,NPARM,MXPARM,CYCMAX,JROUND,ROBUST,prNLL, 1 IFXP,VTP,uVTP,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) Re= PV(NPARM-2) DRMSD= DSE*DSQRT(DFLOAT(NDGF)/DFLOAT(NTP)) IF(prFIT.GT.0) THEN IF(APSE.LE.0) THEN IF(RREF.GT.0.d0) THEN WRITE(6,622) NNAME,q,p,AREF,NL,DRMSD,0, 1 PV(1),PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) ENDIF IF(RREF.LE.0.d0) THEN WRITE(6,624) NNAME,q,p,NL,DRMSD,0,PV(1), 1 PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) ENDIF ELSE IF(RREF.GT.0.d0) WRITE(6,626) NNAME,q,AREF,NS,NL, 1 DRMSD,(j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) IF(RREF.LE.0.d0) WRITE(6,628) NNAME,q,NS,NL,DRMSD, 1 (j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) ENDIF WRITE(6,662) PV(NPbeta+1),PU(NPbeta+1),PS(NPbeta+1) IF(IFXVMIN.LE.0) WRITE(6,660) PV(NPbeta+3), 1 PU(NPbeta+3),PS(NPbeta+3) ENDIF ENDIF c ... then fix Re again, while freeing De & VMIN (as well as the beta's) IF(IFXDe.LE.0) THEN updateCmADJ= -30 DSEB= DSE IFXP(NPARM-2)= 1 NDGF= NTP - NPbeta IF(IFXDe.LE.0) THEN IFXP(NPARM-1)= 0 NDGF= NDGF-1 ENDIF IF(IFXVMIN.LE.0) THEN IFXP(NPARM)= 0 NDGF= NDGF-1 ENDIF CALL NLLSSRR(NTP,NPARM,MXPARM,CYCMAX,JROUND,ROBUST,prNLL, 1 IFXP,VTP,uVTP,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) IF(IFXDe.LE.0) De= PV(NPbeta+2) IF(IFXVMIN.LE.0) VMIN= PV(NPbeta+3) IF((prFIT.GT.0).OR.(DSE.GT.DSEB*1.01)) THEN IF(DSE.GT.DSEB*1.01) WRITE(6,654) DSEB,DSE DRMSD= DSE*DSQRT(DFLOAT(NDGF)/DFLOAT(NTP)) IF(APSE.LE.0) THEN IF(RREF.GT.0.d0) THEN IF((PSEL.EQ.1).OR.(PSEL.EQ.3)) THEN WRITE(6,618) NNAME,q,AREF,NL,DRMSD,0,PV(1), 1 PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) IF(PSEL.EQ.3) WRITE(6,625) AA,BB ENDIF IF(PSEL.EQ.2) THEN WRITE(6,622) NNAME,q,p,AREF,NL,DRMSD,0,PV(1), 1 PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) ENDIF ENDIF IF(RREF.LE.0.d0) THEN IF((PSEL.EQ.1).OR.(PSEL.EQ.3)) THEN WRITE(6,617) NNAME,q,NL,DRMSD,0,PV(1),PU(1), 1 PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) IF(PSEL.EQ.3) WRITE(6,625) AA,BB ENDIF IF(PSEL.EQ.2) THEN WRITE(6,624) NNAME,q,p,NL,DRMSD,0,PV(1), 1 PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=1,NPbeta) ENDIF ENDIF ELSE !! for Pashov Spline exponent IF(RREF.GT.0.d0) WRITE(6,626) NNAME,q,AREF,NS,NL, 1 DRMSD,(j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) IF(RREF.LE.0.d0) WRITE(6,628) NNAME,q,NS,NL,DRMSD, 1 (j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) ENDIF IF(IFXRe.LE.0) 1 WRITE(6,662) PV(NPARM-2),PU(NPARM-2),PS(NPARM-2) IF(IFXDe.LE.0) 1 WRITE(6,630) PV(NPARM-1),PU(NPARM-1),PS(NPARM-1) IF(IFXVMIN.LE.0) 1 WRITE(6,660) PV(NPARM),PU(NPARM),PS(NPARM) ENDIF ENDIF c ... and finally ... fit to all three of VMIN, De and Re IF(IFXDe.LE.0) updateCmADJ= -30 IFXP(NPbeta+1)= IFXRe IFXP(NPbeta+2)= IFXDe IFXP(NPbeta+3)= IFXVMIN PV(NPbeta+1)= Re PV(NPbeta+2)= De PV(NPbeta+3)= VMIN NDGF= NTP- NPbeta IF(IFXVMIN.LE.0) NDGF= NDGF-1 IF(IFXDE.LE.0) NDGF= NDGF-1 IF(IFXRE.LE.0) NDGF= NDGF-1 CALL NLLSSRR(NTP,NPARM,MXPARM,CYCMAX,IROUND,ROBUST,prNLL,IFXP, 1 VTP,uVTP,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) DRMSD= DSE*DSQRT(DFLOAT(NDGF)/DFLOAT(NTP)) IF(APSE.LE.0) THEN IF(RREF.GT.0.d0) THEN IF((PSEL.EQ.1).OR.(PSEL.EQ.3)) THEN WRITE(6,618) NNAME,q,AREF,NL,DRMSD,0,PV(1),PU(1), 1 PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) IF(PSEL.EQ.3) WRITE(6,625) AA,BB WRITE(7,726) NNAME,q,AREF,NL,DRMSD,De,IFXDe,Re, 1 IFXRe,APSE,NL,q,p,AREF,(PV(j),j=1,NPbeta) ENDIF IF(PSEL.EQ.2) THEN WRITE(6,622) NNAME,q,p,AREF,NL,DRMSD,0,PV(1), 1 PU(1),PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) WRITE(7,722) NNAME,q,p,AREF,NL,DRMSD,De,IFXDe,Re, 1 IFXRe,APSE,NL,q,p,AREF,(PV(j),j=1,NPbeta) ENDIF ENDIF IF(RREF.LE.0.d0) THEN IF((PSEL.EQ.1).OR.(PSEL.EQ.3)) THEN WRITE(6,617) NNAME,q,NL,DRMSD,0,PV(1),PU(1),PS(1) 1 ,DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) IF(PSEL.EQ.3) WRITE(6,625) AA,BB WRITE(7,728) NNAME,q,NL,DRMSD,De,IFXDe,Re,IFXRe, 1 APSE,NL,q,p,-AREF,(PV(j),j=1,NPbeta) ENDIF IF(PSEL.EQ.2) THEN WRITE(6,624) NNAME,q,p,NL,DRMSD,0,PV(1),PU(1), 1 PS(1),DSE,(j-1,PV(j),PU(j),PS(j),j=2,NPbeta) WRITE(7,724) NNAME,q,p,NL,DRMSD,De,IFXDe,Re, 1 IFXRe,APSE,NL,q,p,-AREF,(PV(j),j=1,NPbeta) ENDIF ENDIF ELSE !! now print for Pashov Spline Exponent case IF(RREF.GT.0.d0) WRITE(6,626) NNAME,q,AREF,NS,NL,DRMSD, 1 (j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) IF(RREF.LE.0.d0) WRITE(6,628) NNAME,q,NS,NL,DRMSD, 1 (j,yqPSE(j),j,PV(j),PU(j),PS(j),j=1,NPbeta) ENDIF IF(PSEL.EQ.3) BETA0= PV(1) WRITE(6,662) PV(NPbeta+1),PU(NPbeta+1),PS(NPbeta+1) WRITE(6,630) PV(NPbeta+2),PU(NPbeta+2),PS(NPbeta+2) WRITE(6,660) PV(NPbeta+3),PU(NPbeta+3),PS(NPbeta+3) ccc Print [calc.-obs.] IF(prDIFF.gt.0) WRITE(6,730) (RTP(I),VTP(I),YD(I), 1 YD(I)/uVTP(I),I= 1,NTP) 730 FORMAT(1x,38('==')/3x,2(3x,'RTP',7x,'VTP',5x,'[c-o] [c-o]/unc')/ 1 1x,38('--')/(1x,2(0P,f10.5,F10.2,f8.4,1P,D10.2))) ccc IF(IFXRe.LE.0) Re= PV(NPbeta+1) IF(IFXDe.LE.0) De= PV(NPbeta+2) IF(IFXVMIN.LE.0) VMIN= PV(NPbeta+3) ENDIF c======================================================================= c*** For case of a GPEF potential ... c----------------------------------------------------------------------- IF(PSEL.EQ.4) THEN DO i= 1,NTP ccc uVTP(i)= uVTP(i) betay(i)= VTP(i) ENDDO DO NPbeta= NS+1,NL+1 c*** Loop over expansion orders from NS to NL IFXP(NPbeta+1)= IFXVMIN IFXP(NPbeta+2)= IFXRe IFXP(NPbeta+3)= 1 c.... first, do fully linearized fit to get trial expansion coefficients Req= Re**q NPARM= NPbeta IF(IFXVMIN.LE.0) NPARM= NPARM+ 1 cc IF(IFXRe.LE.0) NPARM= NPARM+ 1 DO i= 1, NTP RTPq= RTP(i)**q yq= (RTPq - Req)/(as*RTPq + bs*Req) yPOW= yq DO j=1,NPbeta yPOW= yPOW*yq DYDP(i,j)= yPOW ENDDO IF(IFXVMIN.LE.0) DYDP(i,NPARM)= 1.d0 ENDDO CALL LLSQF(NTP,NPARM,MXDATA,MXPARM,betay,uVTP,DYDP,YD,PV, 1 PU,PS,CM,DSE) IF(IFXVMIN.LE.0) VMIN= PV(NPbeta+1) IF(prFIT.GT.0) THEN WRITE(6,620) NNAME,q,NL,DSE, 1 (' c',j-1,PV(j),PU(j),PS(j),j= 1,NPbeta) IF(IFXVMIN.LE.0) WRITE(6,660) PV(NPbeta+1), 1 PU(NPbeta+1),PS(NPbeta+1) ENDIF c.... then, proceed with fit to non-linear form NPARM= NPbeta+2 PV(NPbeta+2)= Re IFXDe= 1 IF(IFXRe.LE.0) THEN c... If Re is to be free, first optimize it in fit to initial form CALL NLLSSRR(NTP,NPARM,MXPARM,CYCMAX,JROUND,ROBUST, 1 prNLL,IFXP,VTP,uVTP,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) IF(prFIT.GT.0) THEN WRITE(6,636) WRITE(6,658) q,NPbeta,Rref,DSE,(i-1,PV(i),PU(i), 1 PS(i),i=1,NPbeta) IF(IFXVMIN.LE.0) WRITE(6,660) PV(NPbeta+1), 1 PU(NPbeta+1),PS(NPbeta+1) WRITE(6,662) PV(NPbeta+2),PU(NPbeta+2), 1 PS(NPbeta+2) ENDIF ENDIF IF(NPbeta.GE.2) THEN DO j=2, NPbeta PV(j)= PV(j+1)/PV(1) ENDDO ENDIF c ... IFXDe is a flag indicating fit to final c0*y**2(1 + c1*y + ... ) IFXDe= 0 CALL NLLSSRR(NTP,NPARM,MXPARM,CYCMAX,IROUND,ROBUST,prNLL, 1 IFXP,VTP,uVTP,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) WRITE(6,658) q,NPbeta-1,DSE,(i-1,PV(i),PU(i),PS(i), 1 i=1,NPbeta) IF(IFXVMIN.LE.0) WRITE(6,660) PV(NPbeta+1),PU(NPbeta+1), 1 PS(NPbeta+1) IF(IFXRE.LE.0) WRITE(6,662) PV(NPbeta+2),PU(NPbeta+2), 1 PS(NPbeta+2) WRITE(7,659) 2,PV(1),(i+1,PV(1)*pv(I),i=2,nPbeta) ENDDO !! end of loop over polynomial order NPbeta NPARM= NPARM+ 1 ENDIF c====================================================================== c*** Step inward and check whether the repulsive inner potential wall c has an inflection point or turnover. IF(APSE.LE.0) THEN IF(NTP.GE.MXDATA) THEN WRITE(6,674) NTP,MXDATA GOTO 90 ENDIF RH= RTP(1)*1.0d-2 J= NTP+ 1 RR= RTP(1) RTP(J)= RR CALL DYIDPJ(J,MXDATA,NPARM,IFXP,VV,PV,PU,PS) RB= RR + RH RTP(J)= RB CALL DYIDPJ(J,MXDATA,NPARM,IFXP,VB,PV,PU,PS) INFL= 0 DO I= 1,110 RBB= RB VBB= VB RB= RR VB= VV RR= RR - RH RTP(J)= RR IF(RTP(J).LT.0.01d0) EXIT IF(APSE.GT.0) THEN c*** Define SAS value required by dyidpj for this additional point DO m= 1,NPbeta c ... gotta generate xPSE(J) for r= RTP(J) to alow SAS value to be generated XPSE(J)= (RTP(J)**q - AREFq)/(RTP(J)**q + AREFq) SAS(J,m)= Scalc(xPSE(J),m,NPbeta,yqPSE,rKL,MXDATA) ENDDO ENDIF CALL DYIDPJ(J,MXDATA,NPARM,IFXP,VV,PV,PU,PS) cc !!! temporary writes for testing. cc write(6,666) rr,vv, (vv-vb)/RH, (vv - 2.d0*vb + vbb)/RH**2 cc666 format(f8.4,3f16.4) cc !!! end of temporary writes for testing. IF(INFL.GT.0) THEN IF(VV.LT.VB) THEN c... warning of potential inner-wall turnover WRITE(6,672) RB,VB GOTO 90 ENDIF ELSE IF((VV-VB).LE.(VB-VBB)) THEN c... warning of potential inner-wall inflection INFL= 1 WRITE(6,670) RB,VB ENDIF ENDIF ENDDO ENDIF 90 WRITE(6,608) IF((LPPOT.GT.0).AND.(NPR.GT.0)) THEN c** If desired, print fitted potential at NPR points starting at r=RPR1 WRITE(8,800) NNAME, NL, q, p, AREF NPARM= NPbeta+3 VB= 0.d0 SLB= 0.d0 J= NTP+ 1 DO I= 1, NPR RTP(J)= RPR1 + (I-1)*dRPR IF(APSE.GT.0) THEN c*** Define SAS values for this additional point DO m= 1,NPbeta c ... gotta generate xPSE(J) for r= RTP(J) to alow SAS value to be generated XPSE(J)= (RTP(J)**q - AREFq)/(RTP(J)**q + AREFq) SAS(J,m)= Scalc(xPSE(J),m,NPbeta,yqPSE,rKL,MXDATA) ENDDO ENDIF CALL DYIDPJ(J,MXDATA,NPARM,IFXP,VV,PV,PU,PS) SL= (VV-VB)/dRPR dSL= (SL-SLB)/dRPR SLB= SL VB= VV IF(I.EQ.1) WRITE(8,802) RTP(J),VV IF(I.EQ.2) WRITE(8,802) RTP(J),VV, SL IF(I.GT.2) WRITE(8,802) RTP(J),VV, SL,dSL ENDDO ENDIF 800 FORMAT(/2x,A4,' potential for NL=',i3,' q=',I2,' p=',i2, 1 ' Rref=',f8.5/' R V(r)',10x,' dV(r)', 2 8x,'d2V(r)'/2x,29('--') ) 802 FORMAT( f9.4,1pd14.6,2d12.4) GOTO 20 c----------------------------------------------------------------------- 616 FORMAT(//'!!! WARNING !!! Should set p=',i2,' > [{MMLR(',i2, 1 ')=',I2,'} - {MMLR(',i1,')=',i2,'}] CAUTION !!!'/) 619 FORMAT(' Linearized fit uses beta(INF)=',f12.8) 620 FORMAT(/' Linearized ',A4,'{q=',i2,'; NL=',i2,'} fit yields DSE= 1',1Pd9.2/(3x,a4,'_{',i2,'}=',d19.11,' (+/-',d8.1,') PS=', 2 d8.1)) 621 FORMAT(/' Linearized ',A4,'{q=',i2,', p=',i2,'; Rref=',f5.2, 1 '; NL=',i2,'} fit yields DSE=',1Pd9.2/(3x,a4,'_{',i2,'}=', 2 d19.11,' (+/-',d8.1,') PS=',d8.1)) 623 FORMAT(/' Linearized ',A4,'{q=',i2,', p=',i2,'; Rref=Re; NL=', 1 i2,'} fit yields DSE=',1Pd9.2/(3x,a4,'_{',i2,'}=',d19.11, 2 ' (+/-',d8.1,') PS=',d8.1)) 617 FORMAT(/' Direct fit to ',A4,'{q=',i2,'; Rref= Re ; NL=',I2, 1 '} potential:',8x,'dd=',1Pd12.5/ 2 ' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1,3x,'DSE=', 3 D9.2/(' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1)) 618 FORMAT(/' Direct fit to ',A4,'{q=',i2,'; Rref=',f6.3,'; NL=',I2, 1 '} potential:',8x,'dd=',1Pd12.5/ 2 ' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1,3x,'DSE=', 3 D12.5/(' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1)) 622 FORMAT(/' Direct fit to ',A4,'{q=',i2,', p=',i2,'; Rref=',f5.2, 1 '; NL=',I2,'} potential: dd=',1Pd12.5/ 2 ' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1,3x,'DSE=', 3 D12.5/(' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1)) 625 FORMAT(8x,'AA=',1Pd20.12,' BB=',d20.12) 722 FORMAT(/' Direct fit to ',A4,'{q=',i2,', p=',i2,'; Rref=',f5.2, 1 '; NL=',I2,'} potential: dd=',1Pd12.5/d20.12,I3,9x, 2 '% De IFXDe'/d20.12,I3,9x,'% Re IFXRe'//2I3,2I4,D11.2,7x, 3 '% APSE Nbeta nQB nPB RREF'/(d20.12,' 0')) 726 FORMAT(/' Direct fit to ',A4,'{q=',i2,', Rref=',f5.2,'; NL=', 1 I2,'} potential:',7x,'dd=',1Pd12.5/d20.12,I3,9x,'% De IFXDe'/ 1 d20.12,I3,9x,'% Re IFXRe'//2I3,2I4,D11.2,7x,'% APSE Nbeta nQB nP 3B RREF'/(d20.12,' 0')) 624 FORMAT(/' Direct fit to ',A4,'{q=',i2,', p=',I2,'; Rref= Re ; NL= 1',I2,'} potential: dd=',1Pd12.5/ 2 ' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1,3x'DSE=', 3 D9.2/(' beta_{',i2,'}=',d20.12,' (+/-',d8.1,') PS=',d8.1)) 724 FORMAT(/' Direct fit to ',A4,'{q=',i2,', p=',I2,'; Rref= Re ; NL 1=',I2,'} potential: dd=',1Pd12.5/d20.12,I3,9x,'% De IFXDe'/ 2 d20.12,I3,9x,'% Re IFXRe'//2I3,2I4,D11.2,7x,'% APSE Nbeta nQB nP 3B RREF'/(d20.12,' 0')) 728 FORMAT(/' Direct fit to ',A4,'{q=',i2,', Rref= Re ; NL=',I2, 1 '} potential:',7x,'dd=',1Pd12.5/d20.12,I3,9x,'% De IFXDe'/ 2 d20.12,I3,9x,'% Re IFXRe'//2I3,2I4,D11.2,7x,'% APSE Nbeta nQB nP 3B RREF'/(d20.12,' 0')) 626 FORMAT(/' Direct fit to ',A4,'{q=',i2,'; Rref=',f5.2,' ; NS=',i2, 1 ', NL=',I2,'} potential: dd=',1Pd9.2/(' yqPSE{',i2,'}=', 20PF11.7,' beta_{',i2,'}=',1Pd18.10,'(+/-',d8.1,') PS=',d8.1)) 628 FORMAT(/' Direct fit to ',A4,'{q=',i2,'; Rref= Re ; NS=',i2, 1 ', NL=',I2,'} potential: DSE=',1Pd9.2/(' yqPSE{',i2,'}=', 20PF11.7,' beta_{',i2,'}=',1Pd18.10,'(+/-',d8.1,') PS=',d8.1)) 630 FORMAT(8x,'De =',f13.6,' (+/-',f12.6,') PS=',1pd8.1) 632 FORMAT(' Use exponent expansion variable: y_',I1,'(r)= [r^',I1, 1 ' -',f7.4,'^',I1,']/[r^',I1,' +',f7.4,'^',I1,']' ) 633 FORMAT(' Use exponent expansion variable: y_',I1,'(r)= [r^',I1, 1 ' - Re^',I1,']/[r^',I1,' + Re^',I1,']' ) 634 FORMAT(' MLR polynomial exponent function is:'/29x,'beta(R)= betaI 1NF*y_',I1,' + (1-y_',I1,')*Sum{beta_i*[y_',I1,'q]^i}') 636 FORMAT(/' First perform full non-linear GPEF fit without taking ou 1t common factor of c_0') 644 FORMAT(' Update beta_0 from',f11.6,' to',f11.6,' by', 1 1Pd9.1,' : DSE=',1PD8.1) 646 FORMAT(' !!! CAUTION !!! Iteration to optimize beta(0) not conve 1rged after',i3,' tries') 648 FORMAT(' Converge on beta_0=',f11.6,' Next change=',1Pd9.1) 650 FORMAT(' Use Pashov natural spline exponent based on', i4,' yq^{re 1f} values for r < r_e'/41x,'and',i4,' yq^{ref} values for r > r_ 2e') 653 FORMAT(/' Linearized ',A4,'{q=',i2,' p=',i2,'}-APSE treatment yiel 1ds:'/4x,'qpSE',8x,'rPSE',6x,'ypRE'4x,'beta*ypRe beta'/ 2 (F12.6,4f10.6)) 654 FORMAT(/' *** PROBLEM *** freeing De makes DSE increase from', 1 1PD9.2,' to',D9.2) 658 FORMAT(/' Fit to a GPEF{q=',i1,'; N=',I2,'} potential yields:', 1 23x,'DSE=',1Pd9.2/ (5x,'c_{',i2,'} =',d18.10,' (+/-',d8.1,')', 2 4x,'PS=',d8.1)) 659 FORMAT(/(5x,'a_{',i2,'} =',d18.10)) 660 FORMAT(6x,'VMIN =',f13.5,' (+/-',f12.6,') PS=',1pd8.1) 662 FORMAT(8x,'Re =',f13.9,' (+/-',f12.9,') PS=',1pd8.1) 670 FORMAT(1x,39('--')/ ' *** CAUTION *** inner wall has inflection at 1 R=',F6.3,' V=',1PD12.4) 672 FORMAT(28x'and turns over at R=',F6.3,' V=',1PD12.4) 674 FORMAT(/' *** {NTP=',I4,'}.GE.{MXDATA=',i4,'}: array size prevents 1 inner-wall curvature check') 999 STOP END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c======================================================================= SUBROUTINE quadCORR(NCMM,MCMM,NCMMAX,MMLR,De,CmVAL,CmEFF) c======================================================================= c** subroutine to generate and print MLR CmEFF values incorporating c quadratic 'Dattani' corrections to Cm values for both standard 'linear' c and A-F diagonalized uLR(r) functions for MLR potentials c** Return MCMM= NCMM+1 for C9{adj} term for m_1= 3 potentials c======================================================================= INTEGER NCMM,MCMM,NCMMAX,MMLR(NCMMAX) REAL*8 De,CmVAL(NCMMAX),CmEFF(NCMMAX) c----------------------------------------------------------------------- IF(MMLR(1).GT.0) THEN c** For 'normal' inverse-power sum MLR case, with or without damping, c set up Dattani's 'Quadratic-corrected' effective Cm values IF((MMLR(1).EQ.6).AND.(NCMM.GE.4)) THEN c... First, consider C6/C12adj(C14adj) for MMLR(m)={6,8,10,(11),12,14} case IF(MMLR(4).EQ.12) THEN ! explicitly MMLR(4)=12 CmEFF(4)= CmVAL(4)+ 0.25D0*CmVAL(1)**2/De WRITE(6,710) MMLR(4),MMLR(4),CmEFF(4) ENDIF IF(NCMM.GE.5) THEN IF(MMLR(4).EQ.11) THEN ! implicitly MMLR(5)=12 CmEFF(5)= CmVAL(5) + 0.25D0*CmVAL(1)**2/De WRITE(6,710) MMLR(5),MMLR(5),CmEFF(5) IF(NCMM.GE.6) THEN ! implicitly MMLR(6)=14 CmEFF(6)= CmVAL(6)+ 0.5D0*CmVAL(1)*CmVAL(2)/De WRITE(6,710) MMLR(6),MMLR(6),CmEFF(6) ENDIF ENDIF IF(MMLR(4).EQ.12) THEN ! assuming MMLR(5)=14 CmEFF(5)= CmVAL(5) + 0.5D0*CmVAL(1)*CmVAL(2)/De WRITE(6,710) MMLR(5),MMLR(5),CmEFF(5) ENDIF ENDIF ENDIF IF((MMLR(1).EQ.5).AND.(NCMM.GE.4)) THEN c... Then, consider C5/C10adj + C12adj for MMLR(m)={5,6,8,10,12,14} cases CmEFF(4)= CmVAL(4) + 0.25D0*CmVAL(1)**2/De WRITE(6,710) MMLR(4),MMLR(4),CmEFF(4) IF(NCMM.GE.5) THEN ! introduce C12^{adj} CmEFF(5)= CmVAL(5) + 0.25D0*CmVAL(2)**2/De WRITE(6,710) MMLR(5),MMLR(5),CmEFF(5) IF(NCMM.GE.6) THEN ! introduce C14^{adj} CmEFF(6)= CmVAL(6) + 0.5D0*CmVAL(2)*CmVAL(3)/De WRITE(6,710) MMLR(6),MMLR(6),CmEFF(6) ENDIF ENDIF ENDIF IF((MMLR(1).EQ.4).AND.(NCMM.GE.3)) THEN c... Then, consider C4/C8adj + C12adj for MMLR(m)={4,6,8,10,12,14} cases CmEFF(3)= CmVAL(3) + 0.25D0*CmVAL(1)**2/De WRITE(6,712) MMLR(3),MMLR(3),CmEFF(3) IF(NCMM.GE.4) THEN ! implicitly MMLR(4)=10 CmEFF(4)= CmVAL(4) + 0.5D0*CmVAL(1)*CmVAL(2)/De WRITE(6,710) MMLR(4),MMLR(4),CmEFF(4) IF(NCMM.GE.5) THEN ! implicitly MMLR(5)=12 CmEFF(5)= CmVAL(5) + 0.5D0*CmVAL(1)*CmVAL(3)/De 1 + 0.25D0*CmVAL(2)**2/De WRITE(6,710) MMLR(5),MMLR(5),CmEFF(5) IF(NCMM.GE.6) THEN ! implicitly MMLR(6)=14 CmEFF(6)= CmVAL(6)+ 0.5D0*CmVAL(2)*CmVAL(3)/De 1 + 0.5D0*CmVAL(1)*CmVAL(4)/De WRITE(6,710) MMLR(6),MMLR(6),CmEFF(6) ENDIF ENDIF ENDIF ENDIF IF((MMLR(1).EQ.3).AND.(NCMM.GE.2)) THEN c... Then, consider C3/C6adj + C9adj for MMLR(m)={3,6,8,(9),10,12,14} cases CmEFF(2)= CmVAL(2) + 0.25D0*CmVAL(1)**2/De WRITE(6,712) MMLR(2),MMLR(2),CmEFF(2) IF(NCMM.GE.3) THEN ! introduce C9adj & MMLR=9 MCMM= NCMM + 1 MMLR(MCMM)= 9 CmEFF(MCMM)= 0.5d0*CmVAL(1)*CmEFF(2)/De WRITE(6,714) MMLR(MCMM),CmEFF(MCMM) IF(NCMM.GE.5) THEN ! implicitly MMLR(5)=12 CmEFF(5)= CmVAL(5) + 0.5D0*CmVAL(1)*CmEFF(MCMM)/De 1 + 0.25D0*CmEFF(2)**2/De WRITE(6,710) MMLR(5),MMLR(5),CmEFF(5) IF(NCMM.GE.6) THEN ! implicitly MMLR(6)=14 CmEFF(6)= CmVAL(6)+ 0.5D0*CmEFF(2)*CmVAL(3)/De WRITE(6,710) MMLR(6),MMLR(6),CmEFF(6) ENDIF ENDIF ENDIF ENDIF ENDIF c======================================================================= c c** End of CmEFF= Cm + CmADJ setup for non-AF case =================== 710 Format(" 'Quadratic correction' for C",I2,'(MLR) yields', 1 6x,'C',I2,'{adj}=',1PD15.8) 712 Format(" 'Quadratic correction' for C",I1,'(MLR) yields', 1 7x,'C'I1,'{adj}=',1PD15.8) 714 Format(" 'Quadratic corrn' for MLR(m_1=3) introduces C", 1 I1,'(',A4,',adj) =',1PD15.8) 716 Format(" 'Quadratic correction' for C",I1,'(Sigma) yields C', 1 I1,'(Sigma,adj)=',1PD15.8) 718 Format(" 'Quadratic correction' for C",I1,'(^3Pi) yields C', 1 I1,'(^3Pi,adj) =',1PD15.8) 720 Format(" 'Quadratic correction' for C",I1,'(^1Pi) yields C', 1 I1,'(^1Pi,adj) =',1PD15.8) c========================================================================= IF(MMLR(1).LE.0) THEN c** implement Quadratic 'Dattani' MLR corrections for AF cases IF(MMLR(1).GE.-1) THEN !! first for the 2x2 cases ... CmEFF(4)= CmVAL(4) + 0.25*CmVAL(2)**2/De CmEFF(5)= CmVAL(5) + 0.25*CmVAL(3)**2/De WRITE(6,716) MMLR(4),MMLR(4),CmEFF(4) WRITE(6,718) MMLR(5),MMLR(5),CmEFF(5) c* prepare C9{adj} coefficients for addition to chosen root MMLR(8)= 9 !! These terms added just MMLR(9)= 9 !! before exit from AFdiag Cmeff(8)= 0.5*CmVAL(2)*CmEFF(4)/De WRITE(6,714) MMLR(8),'Sigm',CmEFF(8) Cmeff(9)= 0.5*CmVAL(3)*CmEFF(5)/De WRITE(6,714) MMLR(9),'^3Pi',CmEFF(9) ENDIF IF(MMLR(1).LE.-2) THEN !! now for the 3x3 cases ... CmEFF(5)= CmVAL(5) + 0.25*CmVAL(2)**2/De WRITE(6,716) MMLR(5),MMLR(5),CmEFF(5) CmEFF(6)= CmVAL(6) + 0.25*CmVAL(3)**2/De WRITE(6,720) MMLR(6),MMLR(6),CmEFF(6) CmEFF(7)= CmVAL(7) + 0.25*CmVAL(4)**2/De WRITE(6,718) MMLR(7),MMLR(7),CmEFF(7) c* prepare C9{adj} coefficients for addition to chosen root MMLR(11)= 9 !! These terms added just MMLR(12)= 9 !! before exit from AFdiag MMLR(13)= 9 Cmeff(11)= 0.5*CmVAL(2)*CmEFF(5)/De IF(MMLR(1).EQ.-2) WRITE(6,714) MMLR(11),'Sigm',CmEFF(11) Cmeff(12)= 0.5*CmVAL(3)*CmEFF(6)/De IF(MMLR(1).EQ.-3) WRITE(6,714) MMLR(12),'^3Pi',CmEFF(12) Cmeff(13)= 0.5*CmVAL(4)*CmEFF(7)/De IF(MMLR(1).EQ.-4) WRITE(6,714) MMLR(13),'^1Pi',CmEFF(13) ENDIF ENDIF RETURN END c23456789012345678901234567890123456789012345678901234567890123456789012 c*********************************************************************** SUBROUTINE DYIDPJ(IDAT,NDATA,NPTOT,IFXP,YC,PV,PD,PS) c** Subroutine to calculate potential function value YC at distance c RDIST= RTP(IDAT), and its partial derivatives w.r.t. the various c potential parameters. If IDAT.LE.1 generate a new set of c internal potential variables, while if IDAT > 1 use SAVED values c... [Must ensure that calculations based on the current UPDATED PV(j)] c----------------------------------------------------------------------- c********* Introduced AFdiag) 22 December 2015 ********** c********* Last updated 27 January 2015 ********** c----------------------------------------------------------------------- cc INCLUDE 'arrsizes.h' c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c** This 'Block' Data Utility routine with the name 'arrsizes.h', that c governs array dimensioning in program betaFIT MUST be installed c under this name in the same (sub)directory containing the FORTRAN c file(s) for this Program when it is being compiled, or incorporated c into the program wherever dimensioning is required (as it is below) c----------------------------------------------------------------------- INTEGER MXDATA, LMAX, MXPARM, NCMMAX PARAMETER (MXDATA= 1001, LMAX= MXDATA, MXPARM=50, NCMMAX= 25) c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INTEGER i,j,m,IDAT, IFXP(MXPARM),JFXRe,JFXDe,JFXVMIN,NPARM 1 ,NDATA,NPTOT !! 2 param unused here but reqd. by NNLSSRR REAL*8 YC,PV(MXPARM),PD(MXPARM),PS(MXPARM),DM(NCMMAX),DMP(NCMMAX), 1 DMPP(NCMMAX),BETAsm,RTPp,RTPq,Rep,AREF,AREFp,AREFq,ype,dype,dyqe, 2 Scalc,betaINF,yp,yq,yPOW,XP,XPW,DER,TCM,UM,TTMM,DERP,DERPe,DSUM, 3 ypRe,dypRedRe,betaRe,yRePOW,FCT,FCT2,ULR,ULRe,dULRdr, 4 dULRedRe,d2ULRe,DbetaRe,dAAdbi,dAAdRe,dBBdRe,DDER,T0, 5 T0P,T1,RE3,RE6,RE8,RTP3,RTP6,RTP8,RDIST,C6adj,C9adj,BETAN,beta0, 6 C1tst,VLIM REAL*8 DEIGM1(1,1),DEIGM3(1,1),DEIGM5(1,1),DEIGR(1,1),DEIGRe(1,1) c----------------------------------------------------------------------- REAL*8 Re,De,VMIN,RREF,AA,BB,as,bs,rhoAB,BETA, CmVAL(NCMMAX), 1 CmEFF(NCMMAX),dULRdCm(NCMMAX),RTP(MXDATA),SAS(MXDATA,MXPARM) INTEGER PSEL,IFXRe,IFXDe,IFXVMIN,sVSR2,IDSTT,NCMM,MCMM,q,p, 1 NPbeta,APSE,MMLR(NCMMAX),IFXCm(NCMMAX),updateCmADJ COMMON /DATABLK/Re,De,VMIN,RREF,AA,BB,as,bs,rhoAB,CmVAL,CmEFF, 1 dULRdCm,RTP,SAS,BETA,PSEL,IFXRe,IFXDe,IFXVMIN,sVSR2,IDSTT,NCMM, 2 MCMM,q,p,NPbeta,APSE,MMLR,IFXCm,updateCmADJ c----------------------------------------------------------------------- SAVE JFXRe,JFXDe,JFXVMIN, AREF,AREFp,AREFq,Rep,C6adj,C9adj, 1 betaINF,BETA0, ULRe,dULRedRe c======================================================================= IF(ABS(IDAT).LE.1) THEN JFXRe= IFXP(NPbeta+1) JFXDe= IFXP(NPbeta+2) JFXVMIN= IFXP(NPbeta+3) ENDIF RDIST= RTP(IDAT) DO j=1,MXPARM PD(j)= 0.d0 ENDDO c======================================================================= c** For case of an EMO(q) potential c----------------------------------------------------------------------- IF(PSEL.EQ.1) THEN IF(ABS(IDAT).LE.1) THEN IF(JFXRe.LE.0) Re= PV(NPbeta+1) IF(JFXDe.LE.0) De= PV(NPbeta+2) IF(JFXVMIN.LE.0) VMIN= PV(NPbeta+3) AREF= RREF IF(RREF.LE.0) AREF= Re AREFp= AREF**p ENDIF RTPp= RDIST**p yp= (RTPp - AREFp)/(RTPp + AREFp) yPOW= 1.d0 BETAsm= PV(1) DSUM= 0.d0 IF(NPbeta.GE.2) THEN DO j= 2,NPbeta IF(RREF.LE.0.d0) DSUM= DSUM+ (j-1)*PV(j)*yPOW yPOW= yPOW*yp BETAsm= BETAsm+ yPOW*PV(j) ENDDO ENDIF XP= DEXP(-BETAsm*(RDIST- Re)) YC= De*(1.d0 - XP)**2 + VMIN DER= 2.d0*De*(1.d0- XP)*XP DERP= DER*(RDIST- Re) DO j= 1, NPbeta PD(j)= DERP DERP= DERP*yp ENDDO c** If appropriate, also get partial derivative w.r.t. Re IF(JFXRE.LE.0) THEN IF(RREF.LE.0.d0) BETAsm= BETAsm +(RDIST- Re)*DSUM*0.5d0* 1 (p/Re)*(1.d0-yp**2) PD(NPbeta+1)= -DER*BETAsm ENDIF c** If appropriate, also get partial derivative w.r.t. De & VMIN IF(JFXDe.LE.0) PD(NPbeta+2)= (1.d0- XP)**2 IF(JFXVMIN.LE.0) PD(NPbeta+3)= 1.d0 ENDIF c======================================================================= c For the case of an MLR(q,p) potential ... c----------------------------------------------------------------------- IF(PSEL.EQ.2) THEN IF(ABS(IDAT).LE.1) THEN c*** Begin with calculations at r=Re to get uLR(Re) and betaINF IF(JFXRe.LE.0) Re= PV(NPbeta+1) IF(JFXDe.LE.0) De= PV(NPbeta+2) IF(JFXVMIN.LE.0) VMIN= PV(NPbeta+3) AREF= RREF IF(RREF.LE.0.d0) AREF= Re AREFp= AREF**p AREFq= AREF**q Rep= Re**p c!! update Cm{adj} values at beginning of each iteration IF(updateCmADJ.LE.0) CALL quadCORR(NCMM,MCMM,NCMMAX,MMLR, 1 De,CmVAL,CmEFF) updateCmADJ= updateCmADJ+ 1 C!!! IF(MMLR(1).LE.0) THEN c** for Aubert-Frecon 2x2 or 3x3 diagonalization eigenvalues VLIM= VMIN+ De !!!! 0.d0 CALL AFdiag(Re,VLIM,NCMM,NCMMax,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULRe,dULRdCm,dULRedRe) ELSE c** For 'normal' inverse-power sum MLR/MLJ case, with or without damping CALL dampF(Re,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,DM, 1 DMP,DMPP) ULRe= 0.d0 dULRedRe= 0.d0 DO m= 1,MCMM T0= CmEFF(m)/Re**MMLR(m) IF(rhoAB.GT.0.d0) T0= T0*DM(m) ULRe= ULRe + T0 dULRedRe= dULRedRe - MMLR(m)*T0/Re !! corrected!! IF(rhoAB.GT.0.d0) dULRedRe= dULRedRe + 1 T0*DMP(m)/DM(m) ENDDO ENDIF betaINF= DLOG(2.d0*De/ULRe) IF((IDSTT.GT.1).AND.(sVSR2.EQ.-1)) THEN c** For MLR with DS(s=-1/2) damping constrained to have C1= Z1*Z2= IDSTT c then fix the geta(N) value ... BETA0= 0.d0 DO m= 1,MCMM BETA0= BETA0 + CmEFF(m)*DSQRT(3.69d0*rhoAB/ 1 MMLR(m))**(2*MMLR(m)-1) ENDDO C1tst= BETA0/uLre BETAN= DLOG(DSQRT(4.d0*De*IDSTT*116140.97d0)/BETA0) BETA0= DLOG(DSQRT(IDSTT*116140.97d0/De)*uLRe/BETA0) C1tst= De*(DEXP(BETA0)*C1tst)**2 c*** For constrained C1/r case, now define beta(N) BETAN= -0.5D0*BETAN*(-1)**NPbeta T0= 1.d0 DO j= NPbeta-1, 1, -1 BETAN= BETAN + T0*PV(j) T0= -T0 ENDDO PV(NPbeta)= BETAN IFXP(NPbeta)= 1 WRITE(6,666) C1tst/116140.97d0, NPbeta,PV(NPbeta) 666 FORMAT(/' Test C1(sr)=',1Pd10.3,' PV(',i3,')=',1PD16.8/) ENDIF ENDIF c*** END of calculations at r=Re for IDAT=1 to get uLR(Re) and betaINF c--------- now proceed with calcn. for IDAT.ge.1 points -------------- c... First, generate the MLR exponent variables & function ............. RTPp= RDIST**p RTPq= RDIST**q yp= (RTPp - AREFp)/(RTPp + AREFp) yq= (RTPq - AREFq)/(RTPq + AREFq) ype= (RTPp - Rep)/(RTPp + Rep) IF(APSE.GT.0) THEN c*** Case of Pashov natural spline exponent ............................ c... use a spline through the exponent values defined by the input c points to generate values of that exponent at the desired c spline-definition points XP= 0.d0 DO J= 1, NPbeta PD(J)= SAS(IDAT,J) XP= XP + PV(J)*PD(J) ENDDO XP= XP*ype !! hey ... I had dropped the ype here !!! ELSE c... For conventional constrained polynomial exponent function yPOW= 1.d0 - yp BETAsm= PV(1)*yPOW DSUM= 0.d0 IF(NPbeta.GE.2) THEN DO j= 2,NPbeta IF(RREF.LE.0.d0) DSUM= DSUM + PV(j)*(j-1)*yPOW yPOW= yPOW*yq BETAsm= BETAsm+ yPOW*PV(j) ENDDO ENDIF XP= (BETAsm + betaINF*yp)*ype !! total exponent ENDIF c... Now, generate uLR(r) for actual distance ... RDIST=RTP(i) ........ IF(MMLR(1).LE.0) THEN c** For Aubert-Frecon 2x2 or 3x3 diagonalization uLR(r) CALL AFdiag(RDIST,VLIM,NCMM,NCMMax,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULR,dULRdCm,dULRdR) ELSE c** For normal inverse-power sum MLR/MLJ case, with or without damping CALL dampF(RDIST,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,DM, 1 DMP,DMPP) ULR= 0.d0 DO m= 1,MCMM T0= CmEFF(m)/RDIST**MMLR(m) IF(rhoAB.GT.0.d0) T0= T0*DM(m) ULR= ULR + T0 ENDDO ENDIF XPW= DEXP(-XP) * ULR/ULRe YC= De*(1.d0 - XPW)**2 + VMIN DER= 2.d0*De*(1.d0- XPW)*XPW yPOW= DER*ype*(1.d0- yp) IF(APSE.GT.0) THEN c... finalize derivative w.r.t. exponent beta-function spline points ... DO J= 1,NPbeta PD(J)= PD(J)*DER * ype !! I had forgotten ype ENDDO ELSE c... finalize derivative w.r.t. exponent polynomial coefficient .... DO j= 1,NPbeta PD(j)= yPOW yPOW= yPOW*yq ENDDO ENDIF c** If appropriate, also get partial derivative w.r.t. De & VMIN IF(JFXDe.LE.0) PD(NPbeta+2)= (1.d0- XPW)**2 + DER*ype*yp/De IF(JFXVMIN.LE.0) PD(NPbeta+3)= 1.d0 c** If appropriate, also get partial derivative w.r.t. Re IF(JFXRe.LE.0) THEN dype= -0.5d0*(p/RE)*(1.d0 - yp**2) dyqe= -0.5d0*(q/RE)*(1.d0 - yq**2) PD(NPbeta+1)= DER*(dype*XP/ype 1 + (1.d0 - ype*yp)*dULRedRe/ULRe) IF(RREF.LE.0.d0) PD(NPbeta+1)= PD(NPbeta+1) + 1 DER*(dype*(betaINF - BETAsm/(1.d0-yp))+ dyqe*DSUM) ENDIF ENDIF c======================================================================= c** For the case of a DELR_p potential c----------------------------------------------------------------------- IF(PSEL.EQ.3) THEN IF(ABS(IDAT).LE.1) THEN IF(JFXRe.LE.0) Re= PV(NPbeta+1) IF(JFXDe.LE.0) De= PV(NPbeta+2) IF(JFXVMIN.LE.0) VMIN= PV(NPbeta+3) AREF= RREF IF(RREF.LE.0.d0) AREF= Re AREFp= AREF**p ULRe= 0.d0 dULRedRe= 0.d0 d2ULRe= 0.d0 c----------------------------------------------------------------------- c** Evaluate uLR and its first 2 derivs. w,r,t, r, at r=Re ... IF(MMLR(1).LE.0) THEN c** For Aubert-Frecon 2x2 or 3x3 diagonalization uLR(r) CALL AFdiag(Re,VLIM,NCMM,NCMMax,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULRe,dULRdCm,dULRedRe) ELSE CALL dampF(Re,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,DM, 1 DMP,DMPP) ULRe= 0.d0 dULRedRe= 0.d0 d2ULRe= 0.d0 DO m= 1,NCMM T0= CmEFF(m)/Re**MMLR(m) ULRe= ULRe + T0*DM(m) dULRedRe= dULRedRe + T0*(DMP(m)- DM(m)*MMLR(m)/RE) d2ULRe= d2ULRe+T0*(DMPP(m)- 2.d0*MMLR(m)*DMP(m)/RE 1 + MMLR(m)*(MMLR(m)+1.d0)*DM(m)/RE**2) ENDDO Rep= Re**p ypRe= (Rep - AREFp)/(Rep + AREFP) dypRedRe= 2.d0*p*Rep*AREFp/(Re*(Rep + AREFP)**2) IF(RREF.LE.0.d0) dypRedRe= 0.d0 ENDIF c** now get AA & BB their derivs w.r.t. Re & parially w.r.t. beta_i betaRe= PV(1) DbetaRe= 0.d0 !! this is d{beta}/d{y} at r= Re yPOW= 1.d0 IF(NPbeta.GE.1) THEN DO j= 1,NPbeta DbetaRe= DbetaRe + j*PV(j+1)*yPOW yPOW= yPOW*ypRe betaRe= betaRe + PV(J+1)*yPOW !! corrn found by AEM c betaRe= betaRe + yPOW !! error found by AEM ENDDO ENDIF AA= De - ULRe - dULRedRe/betaRe BB= 2.d0*(De - ULRe) - dULRedRe/betaRe dAAdbi= dULRedRe/betaRe**2 dAAdRe= -dULRedRe -d2ULRe/betaRe +dAAdbi*DbetaRe*dypRedRe dBBdRe= dAAdRe - dULRedRe ENDIF c===== end of calcn. for properties at Re performed, for 1'st point ==== RTPp = RDIST**p yp= (RTPp - AREFp)/(RTPp + AREFp) ULR= 0.d0 c ... evaluate uLR at the actual distance ... IF(MMLR(1).LE.0) THEN c** For Aubert-Frecon 2x2 or 3x3 diagonalization uLR(r) CALL AFdiag(RDIST,VLIM,NCMM,NCMMax,MMLR,CmEFF,rhoAB, 1 sVSR2,IDSTT,ULR,dULRdCm,dULRdR) ELSE CALL dampF(RDIST,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,DM, 1 DMP,DMPP) ULR= 0.d0 DO m= 1,NCMM ULR= ULR + DM(m)*CmEFF(m)/RDIST**MMLR(m) ENDDO ENDIF yPOW= 1.d0 yRePOW= 1.d0 BETAsm= PV(1) DSUM= 0.d0 IF(NPbeta.GE.2) THEN DO j= 2,NPbeta IF(RREF.LE.0.d0) DSUM= DSUM+ (j-1)*PV(j)*yPOW yPOW= yPOW*yp BETAsm= BETAsm+ yPOW*PV(j) ENDDO ENDIF XP= DEXP(-BETAsm*(RDIST- Re)) YC= (AA*XP - BB)*XP + De - ULR + VMIN DER= XP*(BB - 2.d0*AA*XP) DERP= DER*(RDIST- Re) DERPe= XP*(XP- 1.d0)*dAAdbi DO j= 1,NPbeta PD(j)= DERP + DERPe DERP= DERP*yp DERPe= DERPe*yp !! another AEM correction ENDDO c** If appropriate, also get partial derivative w.r.t. De & VMIN IF(JFXDe.LE.0) PD(NPbeta+2)= (XP - 2.d0)*XP + 1.d0 IF(JFXVMIN.LE.0) PD(NPbeta+3)= 1.d0 IF(JFXRe.LE.0) THEN c??? For the old sitn. in which r_{ref} \equiv r_e ............... c** If appropriate, also get partial derivative w.r.t. Re PD(NPbeta+1)= BETAsm*(2.d0*XP*AA - BB)*XP 1 + XP*(dAAdRe*XP- dBBdRe) ENDIF ENDIF c======================================================================= c======================================================================= c** For the case of a GPEF(p,as,bs) potential c----------------------------------------------------------------------- IF(PSEL.EQ.4) THEN IF(ABS(IDAT).LE.1) THEN JFXVMIN= IFXP(NPbeta+1) IF(JFXVMIN.LE.0) VMIN= PV(NPbeta+1) JFXRe= IFXP(NPbeta+2) IF(JFXRe.LE.0) Re= PV(NPbeta+2) Rep= Re**p ENDIF RTPp= RDIST**p yp= (RTPp - Rep)/(as*RTPp + bs*Rep) yPOW= yp**2 BETA= PV(1)*yPOW DSUM= 2.d0*PV(1)*yp IF(IFXDe.LE.0) yPOW= PV(1)*yPOW PD(1)= yp**2 IF(NPbeta.GT.1) THEN DO j= 2,NPbeta DSUM= DSUM + (j+1)*PV(j)*yPOW yPOW= yPOW* yp PD(j)= yPOW BETA= BETA+ PV(j)*yPOW ENDDO ENDIF YC= BETA + VMIN IF(JFXVMIN.LE.0) PD(NPbeta+1)= 1.d0 IF(JFXRe.LE.0) THEN PD(NPbeta+2)= -DSUM* (p/Re)*Rep*RTPp*(as+bs)/ 1 (as*RTPp +bs*Rep)**2 ENDIF ENDIF c======================================================================= c%%%%%%% Optional printout for testing cc if(IDAT.eq.1) then cc NPARM= NPbeta cc IF(IFXDe.LE.0) NPARM= NPARM+ 1 cc IF(IFXRe.LE.0) NPARM= NPARM+ 1 cc IF(IFXVMIN.LE.0) NPARM= NPARM+ 1 cc write(8,700) nparm,(PV(i),i=1,nparm) cc write(8,702) (i,i=1,min0(5,nparm)) cc ENDIF cc write(8,704) i,yc,(pd(i),i=1,nparm) cc700 FORMAT(///' Partial derivatives for',i3,' input parameters:'/ cc 1 (1P5D16.8)) cc702 FORMAT(' I YC ',5(' PD(',i1,') ':)) cc704 format(i3,f9.2,1P5D13.5:/(12x,5d13.5:)) c%%%%%%% RETURN END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** c** Asen Pashov's subroutines for constructing spline functions and c their derivatives. double precision function Scalc(x,m,n,XGRID,rKL,LMAX) c** At the position 'x', Scalc is returned as the value of the m'th c of the 'n' Sm(x) function defining a natural cubic spline through the c mesh points located at x= XGRID(x_i), for i=1,n. LMAX specifies the c maximum number of mesh points x= XGRID(x_i) allowed by the calling program c--------------------------------------------------------------------- INTEGER LMAX,I,K,KK,M,N REAL*8 x,y1,y2,XGRID(LMAX),rKL(LMAX,LMAX) k= 0 kk= 0 do i=2,n c... select interval if ((x.gt.XGRID(i-1)).and.(x.le.XGRID(i))) k=i end do if (x.lt.XGRID(1)) then k=2 kk=1 end if if (x.gt.XGRID(n)) then k=n kk=1 end if if(x.eq.XGRID(1)) k=2 y1=XGRID(k-1) y2=XGRID(k) Scalc= 0.d0 IF(kk.eq.0) 1 Scalc= rKL(m,k)*((y1-x)*(((y1-x)/(y1-y2))**2-1)/6)*(y1-y2) 2 + rKL(m,k-1)*((x-y2)*(((x-y2)/(y1-y2))**2-1)/6)*(y1-y2) IF(k.EQ.m) Scalc= Scalc + (y1-x)/(y1-y2) IF(k-1.EQ.m) Scalc= Scalc + (x-y2)/(y1-y2) c... Asen's original coding ... cc Scalc=ndirac(k,m)*A(x,y1,y2)+ndirac(k-1,m)*B(x,y1,y2)+ cc + C(x,y1,y2)*rKL(m,k)+D(x,y1,y2)*rKL(m,k-1) cc else cc Scalc=ndirac(k,m)*A(x,y1,y2)+ndirac(k-1,m)*B(x,y1,y2) cc A=(x1-z)/(x1-x2) cc B=(z-x2)/(x1-x2) cc C=((x1-z)*(((x1-z)/(x1-x2))**2-1)/6)*(x1-x2) cc D=((z-x2)*(((z-x2)/(x1-x2))**2-1)/6)*(x1-x2) c... Asen's original coding ... end c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** double precision function Sprime(x,m,n,XGRID,rKL,LMAX) c** At the position 'x', evaluate the derivative w.r.t. x of the m'th c Sm(x) function contributing the definition of the the natural cubic c spline defined by function values at the n points XGRID(i) [i=1,n] INTEGER i,k,kk,m,n,LMAX REAL*8 x,del,y1,y2,XGRID(LMAX),rKL(LMAX,LMAX) k=0 kk=0 do i=2,n if((x.gt.XGRID(i-1)).and.(x.le.XGRID(i))) k=i enddo if(x.lt.XGRID(1)) then k=2 kk=1 end if if (x.gt.XGRID(n)) then k=n kk=1 end if if (x.eq.XGRID(1)) k=2 y1=XGRID(k-1) y2=XGRID(k) del=y1-y2 Sprime= 0.d0 if(kk.eq.0) Sprime= (del-3.d0*(y1-x)**2/del)*rKL(m,k)/6.d0 + 1 (3.d0*(x-y2)**2/del-del)*rKL(m,k-1)/6.d0 IF(k-1.eq.m) Sprime= Sprime + 1.d0/del IF(k.eq.m) Sprime= Sprime - 1.d0/del ccc if(kk.eq.0) then ccc Sprim=ndirac(k-1,m)/del-ndirac(k,m)/del+ ccc + (del-3*(y1-x)**2/del)*rKL(m,k)/6+ ccc + (3*(x-y2)**2/del-del)*rKL(m,k-1)/6 ccc else ccc Sprim=ndirac(k-1,m)/del-ndirac(k,m)/del ccc end if end c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** subroutine Lkoef(NGRID,XGRID,rKL) c** Call this subroutine with list of the 'NGRID' spline x_i values in c array 'XGRID' with maximum dimension 'LMAX', and it will return the c LMAX x LMAX array of 'rKL' coefficients used for generating the c 'NGRID' S_{NGRID}(x) spline coefficient functions c----------------- Based on nespl subroutine --------------------------- c** CAUTION .. must dimension internal arrays B, INDX & vv @ compilation cc INCLUDE 'arrsizes.h' !! needed only to define LMAX c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c** This 'Block' Data Utility routine with the name 'arrsizes.h', that c governs array dimensioning in program betaFIT MUST be installed c under this name in the same (sub)directory containing the FORTRAN c file(s) for this Program when it is being compiled, or incorporated c into the program wherever dimensioning is required (as it is below) c----------------------------------------------------------------------- INTEGER MXDATA, LMAX, MXPARM, NCMMAX PARAMETER (MXDATA= 1001, LMAX= MXDATA, MXPARM=50, NCMMAX= 25) c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c*** Subroutine by written A. Pashov (2000) ---------------------------- INTEGER I,J,NGRID,INDX(1:LMAX) REAL*8 XGRID(LMAX),rKL(LMAX,LMAX),B(LMAX,LMAX),vv(LMAX), d c ... note vv dimensioned here, but only used in ludcmp !! DO i= 1,LMAX DO j= 1,LMAX rKL(i,j)= 0.d0 B(i,j)= 0.d0 ENDDO ENDDO rKL(1,1)= (XGRID(3)-XGRID(1))/3.d0 rKL(1,2)= (XGRID(3)-XGRID(2))/6.d0 do i= 2,NGRID-3 rKL(i,i-1)= (XGRID(i+1)-XGRID(i))/6.d0 rKL(i,i)= (XGRID(i+2)-XGRID(i))/3.d0 rKL(i,i+1)= (XGRID(i+2)-XGRID(i+1))/6.d0 end do rKL(NGRID-2,NGRID-3)= (XGRID(NGRID-1)-XGRID(NGRID-2))/6.d0 rKL(NGRID-2,NGRID-2)= (XGRID(NGRID)-XGRID(NGRID-2))/3.d0 do i= 1,NGRID-2 B(i,i)= 1.d0/(XGRID(i+1)-XGRID(i)) B(i,i+1)= -1.d0/(XGRID(i+2)-XGRID(i+1))-1.d0/ 1 (XGRID(i+1)-XGRID(i)) B(i,i+2)= 1.d0/(XGRID(i+2)-XGRID(i+1)) end do call ludcmp(rKL,NGRID-2,LMAX,indx,vv,d) do i= 1,NGRID call lubksb(rKL,NGRID-2,LMAX,indx,B(1,i)) end do do i= 1,NGRID-2 do j= 1,NGRID rKL(j,i+1)= B(i,j) end do end do do i= 1,NGRID rKL(i,1)= 0.0d0 rKL(i,NGRID)= 0.0d0 end do end c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE ludcmp(rKL,NGRID,LMAX,indx,vv,d) c** Subroutine taken from "Numerical Recipies", by Press et al. (1986) INTEGER NGRID,LMAX,indx(LMAX),NMAX,i,imax,j,k double precision d,rKL(LMAX,LMAX),vv(LMAX),TINY,aamax,dum,sum PARAMETER (TINY= 1.0e-20) d= 1.d0 do i= 1,NGRID aamax= 0.d0 do j= 1,NGRID if (abs(rKL(i,j)).gt.aamax) aamax= abs(rKL(i,j)) enddo if (aamax.eq.0.) WRITE(6,*) 'singular matrix in ludcmp' vv(i)= 1.d0/aamax enddo do j= 1,NGRID do i= 1,j-1 sum= rKL(i,j) do k= 1,i-1 sum= sum-rKL(i,k)*rKL(k,j) enddo rKL(i,j)= sum enddo aamax= 0.d0 do i= j,NGRID sum= rKL(i,j) do k= 1,j-1 sum= sum-rKL(i,k)*rKL(k,j) enddo rKL(i,j)= sum dum= vv(i)*abs(sum) if (dum.ge.aamax) then imax= i aamax= dum endif enddo if(j.ne.imax)then do k= 1,NGRID dum= rKL(imax,k) rKL(imax,k)= rKL(j,k) rKL(j,k)= dum enddo d= -d vv(imax)= vv(j) endif indx(j)= imax if(rKL(j,j).eq.0.)rKL(j,j)= TINY if(j.ne.NGRID)then dum= 1.d0/rKL(j,j) do i= j+1,NGRID rKL(i,j)= rKL(i,j)*dum enddo endif enddo return END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE lubksb(rKL,NGRID,LMAX,indx,b) c** Subroutine taken from "Numerical Recipies", by Press et al. (1986) INTEGER i,ii,j,ll, NGRID,LMAX,indx(LMAX) double precision rKL(LMAX,LMAX),b(LMAX), sum ii= 0 do i= 1,NGRID ll= indx(i) sum= b(ll) b(ll)= b(i) if (ii.ne.0)then do j= ii,i-1 sum= sum-rKL(i,j)*b(j) enddo else if (sum.ne.0.) then ii= i endif b(i)= sum enddo do i= NGRID,1,-1 sum= b(i) do j= i+1,NGRID sum= sum-rKL(i,j)*b(j) enddo b(i)= sum/rKL(i,i) enddo return END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE NLLSSRR(NDATA,NPTOT,NPMAX,CYCMAX,IROUND,ROBUST,LPRINT, 1 IFXP,YO,YU,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE) c** Program for performing linear or non-linear least-squares fits and c (if desired) automatically using sequential rounding and refitting c to minimize the numbers of parameter digits which must be quoted [see c R.J. Le Roy, J.Mol.Spectrosc. 191, 223-231 (1998)]. 25/03/16 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c COPYRIGHT 1998-2016 by Robert J. Le Roy + c Dept. of Chemistry, Univ. of Waterloo, Waterloo, Ontario, Canada + c This software may not be sold or any other commercial use made + c of it without the express written permission of the author. + c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c** Program uses orthogonal decomposition of the "design" (partial c derivative) matrix for the core locally linear (steepest descent) c step, following a method introduced (to me) by Dr. Michael Dulick. c** If no parameters are free (NPTOT=0), simply return RMS(residuals) as c calculated from the input parameter values {PV(j)}. c** A user MUST SUPPLY subroutine DYIDPJ to generate the predicted c value of each datum and the partial derivatives of each datum w.r.t. c each parameter (see below) from the current trial parameters. c c** On entry: c NDATA is the number of data to be fitted c NPTOT the total number of parameters in the model (.le.NPMAX). c If NPTOT.le.0 , assume YD(i)=YO(i) and calculate the (RMS c dimensionless deviation)=DSE from them & YU(i) c NPMAX is the maximum number of model parameters allowed by current c external array sizes. Should set internal NPINTMX = NPMAX c (may be freely changed by the user). c CYCMAX is the upper bound on the allowed number of iterative cycles c IROUND .ne. 0 causes Sequential Rounding & Refitting to be c performed, with each parameter being rounded at the c |IROUND|'th sig. digit of its local incertainty. c > 0 rounding selects in turn remaining parameter with largest c relative uncertainy c < 0 round parameters sequentially from last to first c = 0 simply stops after full convergence (without rounding). c ROBUST > 0 causes fits to use Watson's ``robust'' weighting c 1/[u^2 +{(c-o)^2}/3]. ROBUST > 1 uses normal 1/u^2 on first c fit cycle and 'robust' on later cycles. c LPRINT specifies the level of printing inside NLLSSRR c if: = 0, no print except for failed convergence. c < 0 only converged, unrounded parameters, PU & PS's c >= 1 print converged parameters, PU & PS's c >= 2 also print parameter change each rounding step c >= 3 also indicate nature of convergence c >= 4 also print convergence tests on each cycle c >= 5 also parameters changes & uncertainties, each cycle c IFXP(j) specifies whether parameter j is to be held fixed c [IFXP > 0] or to be freely varied in the fit [IFXP= 0] c YO(i) are the NDATA 'observed' data to be fitted c YU(i) are the uncertainties in these YO(i) values c PV(j) are initial trial parameter values (for non-linear fits); c should be set at zero for initially undefined parameters. c c** On Exit: c YD(i) is the array of differences [Ycalc(i) - YO(i)] c PV(j) are the final converged parameter values c PU(j) are 95% confidence limit uncertainties in the PV(j)'s c PS(j) are 'parameter sensitivities' for the PV(j)'s, defined such c that the RMS displacement of predicted data due to rounding c off parameter-j by PS(j) is .le. DSE/10*NPTOT c CM(j,k) is the correlation matrix obtained by normalizing variance c /covariance matrix: CM(j,k) = CM(j,k)/SQRT[CM(j,j)*CM(k,k)] c TSTPS = max{|delta[PV(j)]/PS(j)|} is the parameter sensitivity c convergence test: delta[PV(j)] is last change in parameter-j c TSTPU = max{|delta[PV(j)]/PU(j)|} is the parameter uncertainty c convergence test: delta[PV(j)] is last change in parameter-j c DSE is the predicted (dimensionless) standard error of the fit c c NOTE that the squared 95% confidence limit uncertainty in a property c F({PV(j)}) defined in terms of the fitted parameters {PV(j)} (where c the L.H.S. involves [row]*[matrix]*[column] multiplication) is: c [D(F)]^2 = [PU(1)*dF/dPV(1), PU(2)*dF/dPV(2), ...]*[CM(j,k)]* c [PU(2)*dF/dPV(1), PU(2)*dF/dPV(2), ...] c c** Externally dimension: YO, YU and YD .ge. NDATA c PV, PU and PS .ge. NPTOT (say as NPMAX), c CM as a square matrix with column & row length NPMAX c*********************************************************************** INTEGER MXPdim !! internal limit on max # parameters PARAMETER (MXPdim=50) !! must be .GE. external max # NPMAX INTEGER I,J,K,L,IDF,ITER,NITER,CYCMAX,IROUND,ISCAL,JROUND,LPRINT, 1 NDATA,NPTOT,NPMAX,NPARM,NPFIT,JFIX,QUIT,ROBUST, 2 IFXP(NPMAX),JFXP(MXPdim) REAL*8 YO(NDATA), YU(NDATA), YD(NDATA), PV(NPTOT), PU(NPTOT), 1 PS(NPTOT),PSS(MXPdim),PC(MXPdim),PCS(MXPdim),PX(MXPdim), 2 PY(MXPdim),CM(NPMAX,NPMAX), F95(10), 3 RMSR, RMSRB, DSE, TSTPS, TSTPSB, TSTPU, TFACT, S, YC, Zthrd DATA F95/12.7062D0,4.3027D0,3.1824D0,2.7764D0,2.5706D0,2.4469D0, 1 2.3646D0,2.3060D0,2.2622D0,2.2281D0/ IF((NPTOT.GT.NPMAX).OR.(NPTOT.GT.MXPdim).OR.(NPTOT.GT.NDATA) 1 .OR.(NPMAX.GT.MXPdim)) THEN c** If array dimensioning inadequate, print warning & then STOP WRITE(6,602) NPTOT,MXPdim,NPMAX,NDATA STOP ENDIF Zthrd= 0.d0 IF(ROBUST.GE.2) Zthrd= 1.d0/3.d0 TSTPS= 0.d0 RMSR= 0.d0 NITER= 0 QUIT= 0 NPARM= NPTOT DO J= 1, NPTOT PS(J)= 0.d0 JFXP(J)= IFXP(J) IF(IFXP(J).GT.0) NPARM= NPARM- 1 ENDDO NPFIT= NPARM JROUND= IABS(IROUND) c======================================================================= c** Beginning of loop to perform rounding (if desired). NOTE that in c sequential rounding, NPARM is the current (iteratively shrinking) c number of free parameters. 6 IF(NPARM.GT.0) TSTPS= 9.d99 c** TFACT is 95% student t-value for (NDATA-NPARM) degrees of freedom. c [Approximate expression for (NDATA-NPARM).GT.10 accurate to ca. 0.002] TFACT= 0.D0 IF(NDATA.GT.NPARM) THEN IDF= NDATA-NPARM IF(IDF.GT.10) TFACT= 1.960D0*DEXP(1.265D0/DFLOAT(IDF)) IF(IDF.LE.10) TFACT= F95(IDF) ELSE TFACT= 0.D0 ENDIF c====================================================================== c** Begin iterative convergence loop: try for up to 30 cycles DO 50 ITER= 1, CYCMAX ISCAL= 0 NITER= NITER+ 1 DSE= 0.d0 TSTPSB= TSTPS RMSRB= RMSR c** Zero out various arrays 10 IF(NPARM.GT.0) THEN DO I = 1,NPARM c** PSS is the array of Saved Parameter Sensitivities from previous c iteration to be carried into dyidpj subroutine - used in predicting c increment for derivatives by differences. PSS(I)= PS(I) c** PCS is the saved array of parameter changes from previous iteration c to be used (if necessary) to attempt to stablize fit PCS(I)= PC(I) PS(I) = 0.D0 PU(I) = 0.D0 PX(I) = 0.D0 PY(I) = 0.D0 DO J = 1,NPARM CM(I,J) = 0.D0 ENDDO ENDDO ENDIF c c========Beginning of core linear least-squares step==================== c c** Begin by forming the Jacobian Matrix from partial derivative matrix DO I = 1,NDATA c** User-supplied subroutine DYIDPJ uses current (trial) parameter c values {PV} to generate predicted datum # I [y(calc;I)=YC] and its c partial derivatives w.r.t. each of the parameters, returning the c latter in 1-D array PC. See dummy sample version at end of listing. c* NOTE 1: if more convenient, DYIDPJ could prepare the y(calc) values c and derivatives for all data at the same time (when I=1), but only c returned the values here one datum at a time (for I > 1).] c* NOTE 2: the partial derivative array PC returned by DYIDPJ must have c an entry for every parameter in the model, though for parameters c which are held fixed [JFXP(j)=1], those PC(j) values are ignored. CALL DYIDPJ(I,NDATA,NPTOT,JFXP,YC,PV,PC,PSS) IF(NPARM.LT.NPTOT) THEN c** For constrained parameter or sequential rounding, collapse partial c derivative array here DO J= NPTOT,1,-1 IF(JFXP(J).GT.0) THEN c!! First ... move derivative for constrained-parameter POTFIT case cc IF(JFXP(J).GT.1) THEN cc write(6,666) I,J,PC(J),JFXP(J),PC(JFXP(J)) cc PC(JFXP(J))= PC(JFXP(J))+ PC(J) cc666 FORMAT(' For IDAT=',I5,' add PC(',I3,') =',1pD15.8, cc 1 ' to PC(',0pI3,') =',1pD15.8) cc ENDIF c ... now continue collapsing partial derivative array IF(J.LT.NPTOT) THEN DO K= J,NPTOT-1 PC(K)= PC(K+1) ENDDO ENDIF PC(NPTOT)= 0.d0 ENDIF ENDDO ENDIF YD(I)= YC - YO(I) S = 1.D0/YU(I) cc *** For 'Robust' fitting, adjust uncertainties here IF(Zthrd.GT.0.d0) S= 1.d0/DSQRT(YU(I)**2 + Zthrd*YD(I)**2) YC= -YD(I)*S DSE= DSE+ YC*YC IF(NPARM.GT.0) THEN DO J = 1,NPARM PC(J) = PC(J)*S PS(J) = PS(J)+ PC(J)**2 ENDDO CALL QROD(NPARM,NPMAX,NPMAX,CM,PC,PU,YC,PX,PY) ENDIF ENDDO RMSR= DSQRT(DSE/NDATA) IF(NPARM.LE.0) GO TO 60 c c** Compute the inverse of CM CM(1,1) = 1.D0 / CM(1,1) DO I = 2,NPARM L = I - 1 DO J = 1,L S = 0.D0 DO K = J,L S = S + CM(K,I) * CM(J,K) ENDDO CM(J,I) = -S / CM(I,I) ENDDO CM(I,I) = 1.D0 / CM(I,I) ENDDO c c** Solve for parameter changes PC(j) DO I = 1,NPARM J = NPARM - I + 1 PC(J) = 0.D0 DO K = J,NPARM PC(J) = PC(J) + CM(J,K) * PU(K) ENDDO ENDDO c c** Get (upper triangular) "dispersion Matrix" [variance-covarience c matrix without the sigma^2 factor]. DO I = 1,NPARM DO J = I,NPARM YC = 0.D0 DO K = J,NPARM YC = YC + CM(I,K) * CM(J,K) ENDDO CM(I,J) = YC ENDDO ENDDO c** Generate core of Parameter Uncertainties PU(j) and (symmetric) c correlation matrix CM DO J = 1,NPARM PU(J) = DSQRT(CM(J,J)) DO K= J,NPARM CM(J,K)= CM(J,K)/PU(J) ENDDO DO K= 1,J CM(K,J)= CM(K,J)/PU(J) CM(J,K)= CM(K,J) ENDDO ENDDO c c** Generate standard error DSE = sigma^2, and prepare to calculate c Parameter Sensitivities PS IF(NDATA.GT.NPARM) THEN DSE= DSQRT(DSE/(NDATA-NPARM)) ELSE DSE= 0.d0 ENDIF c** Use DSE to get final (95% confid. limit) parameter uncertainties PU c** Calculate 'parameter sensitivities', changes in PV(j) which would c change predictions of input data by an RMS average of DSE*0.1/NPARM YC= DSE*0.1d0/DFLOAT(NPARM) S= DSE*TFACT DO J = 1,NPARM PU(J)= S* PU(J) PS(J)= YC*DSQRT(NDATA/PS(J)) ENDDO c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF((ITER.GT.1).AND.(RMSR.GT.2.0d0*RMSRB).AND.(ISCAL.LE. 3)) 1 THEN c** LeRoy's Marquardt-like attempt to damp changes if RMSR increases ... ISCAL= ISCAL+ 1 IF(LPRINT.GE.0) THEN WRITE(6,620) ITER,RMSR,RMSR/RMSRB,ISCAL 620 FORMAT(' At Iteration',i3,' RMSD=',1PD8.1,' RMSD/RMSDB=',D8.1, 1 " Scale PC by (1/4)**",i1) ccc WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J),J=1,NPTOT) ENDIF DO J= 1,NPTOT PC(J)= 0.25d0*PCS(J) PV(J)= PV(J)- 3.d0*PC(J) ENDDO GOTO 10 ENDIF c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c========End of core linear least-squares step========================== c ... early exit if Rounding cycle finished ... IF(QUIT.GT.0) GO TO 54 c c** Next test for convergence TSTPS= 0.D0 TSTPU= 0.D0 DO J= 1, NPARM TSTPS= MAX(TSTPS,DABS(PC(J)/PS(J))) TSTPU= MAX(TSTPU,DABS(PC(J)/PU(J))) ENDDO IF(LPRINT.GE.4) WRITE(6,604) ITER,RMSR,TSTPS,TSTPU c** Now ... update parameters (careful about rounding) DO J= 1,NPTOT IF(JFXP(J).GT.0) THEN cc IF(JFXP(J).GT.1) THEN !! a special PotFit option c** If this parameter constrained to equal some earlier parameter .... cc PV(J)= PV(JFXP(J)) cc WRITE(6,668) J,JFXP(J),PV(J),ITER cc ENDIF cc668 FORMAT(' Constrain PV('i3,') = PV(',I3,') =',1pd15.8, cc 1 ' on cycle',i3) c** If parameter held fixed (by input or rounding process), shift values c of change, sensitivity & uncertainty to correct label. IF(J.LT.NPTOT) THEN DO I= NPTOT,J+1,-1 PC(I)= PC(I-1) PS(I)= PS(I-1) PU(I)= PU(I-1) ENDDO ENDIF PC(J)= 0.d0 PS(J)= 0.d0 PU(J)= 0.d0 ELSE PV(J)= PV(J)+ PC(J) ENDIF ENDDO IF(LPRINT.GE.5) WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J), 1 J=1,NPTOT) IF(ITER.GT.1) THEN c** New Convergence test requires RMSD to be constant to 1 part in 10^7 c in adjacent cycles (unlikely to occur by accident!) c** Replaces less severe requirement that TSTPS < 1.0 IF(ABS((RMSR/RMSRB)-1.d0).LT.1.d-07) THEN IF(LPRINT.GE.3) WRITE(6,607) ITER, 1 ABS(RMSR/RMSRB-1.d0),TSTPS GO TO 54 ENDIF ENDIF cc CALL FLUSH(6) IF(ROBUST.GT.0) Zthrd= 1.d0/3.d0 50 CONTINUE WRITE(6,610) NPARM,NDATA,ITER,RMSR,TSTPS,TSTPU c** End of iterative convergence loop for (in general) non-linear case. c====================================================================== c 54 IF(NPARM.LT.NPTOT) THEN c** If necessary, redistribute correlation matrix elements to full c NPTOT-element correlation matrix DO J= 1,NPTOT IF(JFXP(J).GT.0) THEN c* If parameter J was held fixed IF(J.LT.NPTOT) THEN c ... then move every lower CM element down one row: DO I= NPTOT,J+1,-1 c ... For K < J, just shift down or over to the right IF(J.GT.1) THEN DO K= 1,J-1 CM(I,K)= CM(I-1,K) CM(K,I)= CM(I,K) ENDDO ENDIF c ... while for K > J also shift elements one column to the right DO K= NPTOT,J+1,-1 CM(I,K)= CM(I-1,K-1) ENDDO ENDDO ENDIF c ... and finally, insert appropriate row/column of zeros .... DO I= 1,NPTOT CM(I,J)= 0.d0 CM(J,I)= 0.d0 ENDDO CM(J,J)= 1.d0 ENDIF ENDDO ENDIF IF(QUIT.GT.0) GOTO 60 IF(NPARM.EQ.NPFIT) THEN c** If desired, print unrounded parameters and fit properties IF(LPRINT.NE.0) THEN WRITE(6,616) NDATA,NPARM,RMSR,TSTPS WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J),J=1,NPTOT) ENDIF ENDIF IF(IROUND.EQ.0) RETURN c** Automated 'Sequential Rounding and Refitting' section: round c selected parameter, fix it, and return (above) to repeat fit. IF(IROUND.LT.0) THEN c ... if IROUND < 0, sequentially round off 'last' remaining parameter DO J= 1, NPTOT IF(JFXP(J).LE.0) THEN JFIX= J ENDIF ENDDO ELSE c ... if IROUND > 0, sequentially round off remaining parameter with c largest relative uncertainty. c ... First, select parameter JFIX with the largest relative uncertainty JFIX= NPTOT K= 0 TSTPS= 0.d0 DO J= 1,NPTOT IF(JFXP(J).LE.0) THEN K= K+1 TSTPSB= DABS(PU(J)/PV(J)) IF(TSTPSB.GT.TSTPS) THEN JFIX= J TSTPS= TSTPSB ENDIF ENDIF ENDDO ENDIF YC= PV(JFIX) CALL ROUND(JROUND,NPMAX,NPTOT,NPTOT,JFIX,PV,PU,PS,CM) JFXP(JFIX)= 1 IF(LPRINT.GE.2) 1 WRITE(6,614) JFIX,YC,PU(JFIX),PS(JFIX),JFIX,PV(JFIX),RMSR NPARM= NPARM-1 IF(NPARM.EQ.0) THEN c** After rounding complete, make one more pass with all non-fixed c parameters set free to get full correct final correlation matrix, c uncertainties & sensitivities. Don't update parameters on this pass! NPARM= NPFIT QUIT= 1 DO J= 1,NPTOT JFXP(J)= IFXP(J) ENDDO c ... reinitialize for derivative-by-differences calculation RMSR= 0.d0 ENDIF GO TO 6 c c** If no parameters varied or sequential rounding completed - simply c calculate DSE from RMS residuals and return. 60 DSE= 0.d0 IF(NDATA.GT.NPFIT) THEN DSE= RMSR*DSQRT(DFLOAT(NDATA)/DFLOAT(NDATA-NPFIT)) ELSE DSE= 0.d0 ENDIF IF(NPFIT.GT.0) THEN IF(LPRINT.GT.0) THEN c** Print final rounded parameters with original Uncert. & Sensitivities IF(QUIT.LT.1) WRITE(6,616) NDATA, NPFIT, RMSR, TSTPS IF(QUIT.EQ.1) WRITE(6,616) NDATA, NPFIT, RMSR DO J= 1, NPTOT IF(JFXP(J).GT.0) THEN c** If parameter held fixed (by rounding process), shift values of c change, sensitivity & uncertainty to correct absolute number label. DO I= NPTOT,J+1,-1 PC(I)= PC(I-1) PS(I)= PS(I-1) PU(I)= PU(I-1) ENDDO PC(J)= 0.d0 PS(J)= 0.d0 PU(J)= 0.d0 ENDIF ENDDO WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J),J=1,NPTOT) ENDIF ENDIF RETURN c 602 FORMAT(/' *** NLLSSRR problem: [NPTOT=',i4,'] > min{MXPdim=', 1 i4,' NPMAX=',i4,', NDATA=',i6,'}') 604 FORMAT(' After Cycle #',i2,': DRMSD=',1PD14.7,' test(PS)=', 1 1PD8.1,' test(PU)=',D8.1) 606 FORMAT(/' Effective',i3,'-cycle Cgce: MAX{|change/unc.|}=', 1 1PD8.1,' < 0.01 DRMSD=',D10.3) 607 FORMAT(/' Full',i3,'-cycle convergence: {ABS(RMSR/RMSRB)-1}=', 1 1PD9.2,' TSTPS=',D8.1) 610 FORMAT(/ ' !! CAUTION !! fit of',i5,' parameters to',I6,' data not 1 converged after',i3,' Cycles'/5x,'DRMS(deviations)=',1PD10.3, 2 ' test(PS) =',D9.2,' test(PU) =',D9.2/1x,31('**')) 612 FORMAT((3x,'PV(',i4,') =',1PD22.14,' (+/-',D8.1,') PS=',d8.1, 1 ' PC=',d9.1)) 614 FORMAT(' =',39('==')/' Round Off PV(',i4,')=',1PD21.13,' (+/-', 1 D9.2,') PS=',d9.2/4x,'fix PV(',I4,') as ',D19.11, 2 ' & refit: DRMS(deviations)=',D12.5) 616 FORMAT(/i6,' data fit to',i5,' param. yields DRMS(devn)=', 1 1PD14.7:' tst(PS)=',D8.1) END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE QROD(N,NR,NC,A,R,F,B,GC,GS) C** Performs ORTHOGONAL DECOMPOSITION OF THE LINEAR LEAST-SQUARES C EQUATION J * X = F TO A * X = B(TRANSPOSE) * F WHERE C J IS THE JACOBIAN IN WHICH THE FIRST N ROWS AND COLUMNS C ARE TRANSFORMED TO THE UPPER TRIANGULAR MATRIX A C (J = B * A), X IS THE INDEPENDENT VARIABLE VECTOR, AND C F IS THE DEPENDENT VARIABLE VECTOR. THE TRANSFORMATION C IS APPLIED TO ONE ROW OF THE JACOBIAN MATRIX AT A TIME. C PARAMETERS : C N - (INTEGER) DIMENSION OF A TO BE TRANSFORMED. C NR - (INTEGER) ROW DIMENSION OF A DECLARED IN CALLING PROGRAM. C NC - (INTEGER) Column DIMENSION OF F DECLARED IN CALLING PROGRAM. C A - (REAL*8 ARRAY OF DIMENSIONS .GE. N*N) UPPER TRIANGULAR C TRANSFORMATION MATRIX. C R - (REAL*8 LINEAR ARRAY OF DIMENSION .GE. N) ROW OF C JACOBIAN TO BE ADDED. C F - (REAL*8 LINEAR ARRAY .GE. TO THE ROW DIMENSION OF THE C JACOBIAN) TRANSFORMED DEPENDENT VARIABLE MATRIX. C B - (REAL*8) VALUE OF F THAT CORRESPONDS TO THE ADDED C JACOBIAN ROW. C GC - (REAL*8 LINEAR ARRAY .GE. N) GIVENS COSINE TRANSFORMATIONS. C GS - (REAL*8 LINEAR ARRAY .GE. N) GIVENS SINE TRANSFORMATIONS. C-------------------------------------------------------------------- C AUTHOR : MICHAEL DULICK, Department of Chemistry, C UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO N2L 3G1 C-------------------------------------------------------------------- INTEGER I,J,K,N,NC,NR REAL*8 A(NR,NC), R(N), F(NR), GC(N), GS(N), B, Z(2) DO 10 I = 1,N Z(1) = R(I) J = I - 1 DO K = 1,J Z(2) = GC(K) * A(K,I) + GS(K) * Z(1) Z(1) = GC(K) * Z(1) - GS(K) * A(K,I) A(K,I) = Z(2) ENDDO GC(I) = 1.D0 GS(I) = 0.D0 IF(DABS(Z(1)).LE.0.D0) GOTO 10 IF(DABS(A(I,I)) .LT. DABS(Z(1))) THEN Z(2) = A(I,I) / Z(1) GS(I) = 1.D0 / DSQRT(1.D0 + Z(2) * Z(2)) GC(I) = Z(2) * GS(I) ELSE Z(2) = Z(1) / A(I,I) GC(I) = 1.D0 / DSQRT(1.D0 + Z(2) * Z(2)) GS(I) = Z(2) * GC(I) ENDIF A(I,I) = GC(I) * A(I,I) + GS(I) * Z(1) Z(2) = GC(I) * F(I) + GS(I) * B B = GC(I) * B - GS(I) * F(I) F(I) = Z(2) 10 CONTINUE RETURN END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE ROUND(IROUND,NPMAX,NPARM,NPTOT,IPAR,PV,PU,PS,CM) c** Subroutine to round off parameter # IPAR with value PV(IPAR) at the c |IROUND|'th significant digit of: [its uncertainty PU(IPAR)] . c** On return, the rounded value replaced the initial value PV(IPAR). c** Then ... use the correlation matrix CM and the uncertainties PU(I) c in the other (NPTOT-1) [or (NPARM-1) free] parameters to calculate c the optimum compensating changes PV(I) in their values. c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c COPYRIGHT 1998 by Robert J. Le Roy + c Dept. of Chemistry, Univ. of Waterloo, Waterloo, Ontario, Canada + c This software may not be sold or any other commercial use made + c of it without the express written permission of the author. + c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ INTEGER IROUND,NPMAX,NPARM,NPTOT,IPAR,I,IRND,KRND REAL*8 PU(NPMAX),PS(NPMAX),PV(NPMAX),CM(NPMAX,NPMAX),CNST, 1 CRND,XRND,FCT,Z0 DATA Z0/0.d0/ CNST= PV(IPAR) XRND= DLOG10(PU(IPAR)) c** If appropriate, base last rounding step on sensitivity (not uncert.) IF((NPARM.EQ.1).AND.(PS(IPAR).LT.PU(IPAR))) XRND= DLOG10(PS(IPAR)) c** First ... fiddle with log's to perform the rounding IRND= INT(XRND) IF(XRND.GT.0) IRND=IRND+1 IRND= IRND- IROUND FCT= 10.D0**IRND CRND= PV(IPAR)/FCT XRND= Z0 c ... if rounding goes past REAL*8 precision, retain unrounded constant IF(DABS(CRND).GE.1.D+16) THEN WRITE(6,601) IROUND,IPAR RETURN ENDIF IF(DABS(CRND).GE.1.D+8) THEN c ... to avoid problems from overflow of I*4 integers ... KRND= NINT(CRND/1.D+8) XRND= KRND*1.D+8 CRND= CRND-XRND XRND= XRND*FCT END IF IRND= NINT(CRND) CNST= IRND*FCT+ XRND c???????????????? c** Zero parameters more aggressively ... if unc. > 2* value if(dabs(PU(IPAR)/PV(IPAR)).GT.2.d0) then cnst= 0.d0 endif c???????????????? c** Now ... combine rounding change in parameter # IPAR, together with c correlation matrix CM and parameter uncertainties PU to predict c changes in other parameters to optimally compensate for rounding off c of parameter-IPAR. Method pointed out by Mary Thompson (Dept. of c Statistics, UW), IF(IPAR.GT.1) THEN XRND= (CNST-PV(IPAR))/PU(IPAR) DO I= 1,NPTOT IF(I.NE.IPAR) THEN PV(I)= PV(I)+ CM(IPAR,I)*PU(I)*XRND ENDIF ENDDO ENDIF PV(IPAR)= CNST RETURN 601 FORMAT(' =',39('==')/' Caution:',i3,'-digit rounding of parameter- 1',i2,' would exceed (assumed) REAL*8'/' ******** precision overf 2low at 1.D+16, so keep unrounded constant') END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** c SUBROUTINE DYIDPJ(I,NDATA,NPTOT,IFXP,YC,PV,PD.PS) c** Illustrative dummy version of DYIDPJ for the case of a fit to a c power series of order (NPTOT-1) in X(i). *** For datum number-i, c calculate and return PD(j)=[partial derivatives of datum-i] w.r.t. c each of the free polynomial coefficients varied in the fit c (for j=1 to NPTOT). ** Elements of the integer array IFXP indicate c whether parameter j is being held fixed [IFXP(j) > 0] or varied in c the fit [IFXP(j).le.0]. If the former, the partial derivative c for parameter j should be PD(j)= 0.0. c===================================================================== c** Use COMMON block(s) to bring in values of the independent variable c [here XX(i)] and any other parameters or variables needeed to c calculate YC and the partial derivatives. c===================================================================== c INTEGER I,J,NDATA,NPTOT,MXDATA,IFXP(NPTOT) c PARAMETER (MXDATA= 501) c REAL*8 RMSR,YC,PV(NPTOT),PD(NPTOT),PS(NPTOT),POWER,XX(MXDATA) c COMMON /DATABLK/XX c===================================================================== c** NOTE BENE(!!) for non-linear fits, need to be sure that the c calculations of YC and PD(j) are based on the current UPDATED PV(j) c values. If other (than PV) parameter labels are used internally c in the calculations, UPDATE them whenever (say) I = 1 . c===================================================================== c POWER= 1.D0 c YC= PV(1) c PD(1)= POWER c DO 10 J= 2,NPTOT c POWER= POWER*XX(I) c YC= YC+ PV(J)*POWER c PD(J)= POWER c 10 CONTINUE c RETURN c END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE LLSQF(NDATA,NPARM,MXDATA,MXPARM,YO,YU,DYDP,YD,PV,PU,PS, 1 CM,DSE) c** Program for performing linear least squares fits using orthogonal c decomposition of the Design (partial derivative) matrix. c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c COPYRIGHT 2000 by Robert J. Le Roy + c Dept. of Chemistry, Univ. of Waterloo, Waterloo, Ontario, Canada + c This software may not be sold or any other commercial use made + c of it without the express written permission of the author. + c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c** This version of the program is designed for the data sets of modest c size where it is convenient to generate and store the complete c partial derivative matrix prior to calling LLSQF. If this is not the c case, subroutine version LLSQFVL, which generates this partial c derivative array one row at a time through calls to a user-supplied c subroutine, should be used. c c** On entry: NDATA is the number of data to be fitted (.le.MXDATA) c NPARM the number of parameters to be varied (.le.MXPARM) c If NPARM.le.0 , assume YD(i)=YO(i) and calculate the c (RMS dimensionless deviation)=DSE from them & YU(i) c MXDATA & MXPARM are array dimension parameters (see below) c Internal array sizes currently assume MXPARM .le. 60 c YO(i) are the NDATA 'observed' data; for iterative c non-linear fits these are: [Y(obs,i) - Y(trial,i)] c YU(i) are the uncertainties in these YO(i) values c DYDP(i,j) is the partial derivative array dYO(i)/dPV(j) c c** On Exit: PV(j) are the fitted parameter values; for iterative c non-linear fits these are the parameter changes c PU(j) are 95% confidence limit uncertainties in the PV(j)'s c PS(j) are 'parameter sensitivities' for the PV(j)'s, defined c such that the RMS displacement of predicted data due to c rounding off parameter-j by PS(j) is .le. DSE/10*NPARM c DSE is predicted (dimensionless) standard error of the fit c YD(i) is the array of differences [YO(i) - Ycalc(i)] c CM(j,k) is the correlation matrix obtained by normalizing c variance/covariance matrix: CM(j,k) = CM(j,k)/SQRT[CM(j,j)*CM(k,k)] c** The squared 95% confidence limit uncertainty in a property F({PV(j)}) c defined in terms of the fitted parameters {PV(j)} is (where the c L.H.S. involves [row]*[matrix]*[column] multiplication): c [D(F)]^2 = [PU(1)*dF/dPV(1), PU(2)*dF/dPV(2), ...]*[CM(j,k)]* c [PU(2)*dF/dPV(1), PU(2)*dF/dPV(2), ...] c** Externally dimension: YO, YU and YD .ge. NDATA (say as MXDATA), c PV, PU and PS .ge. NPARM (say as MXPARM), c DYDP with column length MXDATA and row length .ge. NPARM c CM as a square matrix with column & row length MXPARM c Authors: Robert J. Le Roy & Michael Dulick, Department of Chemistry c U. of Waterloo, Waterloo, Ontario N2L 3G1. Version of: 07/10/00 c*********************************************************************** INTEGER I,J,K,L,IDF,NDATA,MXDATA,NPARM,MXPARM REAL*8 YO(NDATA), YU(NDATA), YD(NDATA), PV(NPARM), PU(NPARM), 1 PS(NPARM), DYDP(MXDATA,NPARM), CM(MXPARM,MXPARM), DSE, 2 PX(60), F95(10), TFACT, S, U DATA F95/12.7062D0,4.3027D0,3.1824D0,2.7764D0,2.5706D0,2.4469D0, 1 2.3646D0,2.3060D0,2.2622D0,2.2281D0/ c IF((NDATA.GT.MXDATA).OR.(NPARM.GT.MXPARM).OR.(NPARM.GT.60) 1 .OR.(NPARM.GT.NDATA)) THEN c** If array dimensioning inadequate, print warning & then STOP WRITE(6,601) NDATA,MXDATA,NPARM,MXPARM STOP ENDIF IF(NPARM.LE.0) THEN c** If no parameters varied - simply calculate RMS deviation = DSE DSE= 0.D0 DO I= 1,NDATA YD(I)= YO(I) DSE= DSE+ (YD(I)/YU(I))**2 ENDDO DSE= DSQRT(DSE/DFLOAT(NDATA)) RETURN ENDIF c** TFACT is 95% student t-value for (NDATA-NPARM) degrees of freedom. c [Approximate expression for (NDATA-NPARM).GT.10 accurate to ca. 0.002] TFACT= 0.D0 IF(NDATA.GT.NPARM) THEN IDF= NDATA-NPARM IF(IDF.GT.10) TFACT= 1.960D0*DEXP(1.265D0/DFLOAT(IDF)) IF(IDF.LE.10) TFACT= F95(IDF) ELSE TFACT= 0.D0 ENDIF DO I = 1,NPARM PS(I) = 0.D0 PU(I) = 0.D0 PX(I) = 0.D0 DO J = 1,NPARM CM(I,J) = 0.D0 ENDDO ENDDO c c** Begin by forming the Jacobian Matrix from the input partial c derivative matrix DYDP. For VERY large data sets, these partial c derivatives may be generated inside this loop (see version LLSQFVL). DO I = 1,NDATA S = 1.D0 / YU(I) U = YO(I) * S DO J = 1,NPARM PV(J) = DYDP(I,J) * S ENDDO CALL QQROD(NPARM,MXPARM,MXPARM,CM,PV,PX,U,PS,PU) ENDDO c c** Compute the inverse of CM CM(1,1) = 1.D0 / CM(1,1) DO I = 2,NPARM L = I - 1 DO J = 1,L S = 0.D0 DO K = J,L S = S + CM(K,I) * CM(J,K) ENDDO CM(J,I) = -S / CM(I,I) ENDDO CM(I,I) = 1.D0 / CM(I,I) ENDDO c c** Solve for parameter values PV(j) DO I = 1,NPARM J = NPARM - I + 1 PV(J) = 0.D0 DO K = J,NPARM PV(J) = PV(J) + CM(J,K) * PX(K) ENDDO ENDDO c c** Get (upper triangular) "dispersion Matrix" [variance-covarience c matrix without the sigma^2 factor]. DO I = 1,NPARM DO J = I,NPARM U = 0.D0 DO K = J,NPARM U = U + CM(I,K) * CM(J,K) ENDDO CM(I,J) = U ENDDO ENDDO c** Generate core of Parameter Uncertainties PU(j) and (symmetric) c correlation matrix CM DO J = 1,NPARM PU(J) = DSQRT(CM(J,J)) DO K= J,NPARM CM(J,K)= CM(J,K)/PU(J) ENDDO DO K= 1,J CM(K,J)= CM(K,J)/PU(J) CM(J,K)= CM(K,J) ENDDO PX(J)= 0.d0 ENDDO c c** Generate differences: YD(i) = [YO(i) - Ycalc(i)] , standard error c DSE = sigma^2, and prepare to calculate Parameter Sensitivities PS DSE= 0.D0 DO I = 1,NDATA S = 1.D0 / YU(I) U = 0.D0 DO J = 1,NPARM PX(J)= PX(J)+ (DYDP(I,J)*S)**2 U = U + DYDP(I,J) * PV(J) ENDDO YD(I) = YO(I) - U DSE= DSE+ (S*YD(I))**2 ENDDO IF(NDATA.GT.NPARM) THEN DSE= DSQRT(DSE/(NDATA-NPARM)) ELSE DSE= 0.d0 ENDIF c** Use DSE to get final (95% confid. limit) parameter uncertainties PU c** Calculate 'parameter sensitivities', changes in PV(j) which would c change predictions of input data by an RMS average of DSE*0.1/NPARM U= DSE*0.1d0/DFLOAT(NPARM) S= DSE*TFACT DO J = 1,NPARM PU(J)= S* PU(J) PS(J)= U*DSQRT(NDATA/PX(J)) ENDDO c RETURN 601 FORMAT(/' *** Dimensioning problems in LLSQF *** (NDATA, MXDATA, N 1PARM, MXPARM) = (',I5,4(' ,',I5),' )') END c*********************************************************************** SUBROUTINE QQROD(N,NR,NC,A,R,F,B,GC,GS) C** Performs ORTHOGONAL DECOMPOSITION OF THE LINEAR LEAST-SQUARES C EQUATION J * X = F TO A * X = B(TRANSPOSE) * F WHERE C J IS THE JACOBIAN IN WHICH THE FIRST N ROWS AND COLUMNS C ARE TRANSFORMED TO THE UPPER TRIANGULAR MATRIX A C (J = B * A), X IS THE INDEPENDENT VARIABLE VECTOR, AND C F IS THE DEPENDENT VARIABLE VECTOR. THE TRANSFORMATION C IS APPLIED TO ONE ROW OF THE JACOBIAN MATRIX AT A TIME. C PARAMETERS : C N - (INTEGER) DIMENSION OF A TO BE TRANSFORMED. C NR - (INTEGER) ROW DIMENSION OF A DECLARED IN CALLING PROGRAM. C NC - (INTEGER) Column DIMENSION OF F DECLARED IN CALLING PROGRAM. C A - (REAL*8 ARRAY OF DIMENSIONS .GE. N*N) UPPER TRIANGULAR C TRANSFORMATION MATRIX. C R - (REAL*8 LINEAR ARRAY OF DIMENSION .GE. N) ROW OF C JACOBIAN TO BE ADDED. C F - (REAL*8 LINEAR ARRAY .GE. TO THE ROW DIMENSION OF THE C JACOBIAN) TRANSFORMED DEPENDENT VARIABLE MATRIX. C B - (REAL*8) VALUE OF F THAT CORRESPONDS TO THE ADDED C JACOBIAN ROW. C GC - (REAL*8 LINEAR ARRAY .GE. N) GIVENS COSINE TRANSFORMATIONS. C GS - (REAL*8 LINEAR ARRAY .GE. N) GIVENS SINE TRANSFORMATIONS. C-------------------------------------------------------------------- C AUTHOR : MICHAEL DULICK, Department of Chemistry, C UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO N2L 3G1 C-------------------------------------------------------------------- INTEGER I,J,K,N,NC,NR REAL*8 A(NR,NC), R(N), F(NR), GC(N), GS(N), B, Z(2) DO I = 1,N Z(1) = R(I) J = I - 1 DO K = 1,J Z(2) = GC(K) * A(K,I) + GS(K) * Z(1) Z(1) = GC(K) * Z(1) - GS(K) * A(K,I) A(K,I) = Z(2) ENDDO GC(I) = 1.D0 GS(I) = 0.D0 IF((Z(1) .GT. 0.d0) .OR. (Z(1) .LT. 0.d0)) THEN IF(DABS(A(I,I)) .LT. DABS(Z(1))) THEN Z(2) = A(I,I) / Z(1) GS(I) = 1.D0 / DSQRT(1.D0 + Z(2) * Z(2)) GC(I) = Z(2) * GS(I) ELSE Z(2) = Z(1) / A(I,I) GC(I) = 1.D0 / DSQRT(1.D0 + Z(2) * Z(2)) GS(I) = Z(2) * GC(I) ENDIF A(I,I) = GC(I) * A(I,I) + GS(I) * Z(1) Z(2) = GC(I) * F(I) + GS(I) * B B = GC(I) * B - GS(I) * F(I) F(I) = Z(2) ENDIF ENDDO RETURN END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE AFdiag(RDIST,VLIM,NCMM,NCMMax,MMLR,Cm,rhoAB,sVSR2, 1 IDSTT,ULR,dULRdCm,dULRdR) c*********************************************************************** c** Aubert-Frecon Potential Model for u_{LR}(r) c*********************************************************************** c** Subroutine to generate, at the onee distance RDIST, an eigenvalue c of the 2x2 or 3x3 long-range interaction matrix described by Eqs.1 c and 10, resp., of J.Mol.Spec.188, 182 (1998) (Aubert-Frecon et al) c** and its derivatives w.r.t. the C_m long-range parameters. c*********************************************************************** c==> Input: r= RDIST, VLIM, NCMM, m=MMLR & Cm's, rhoAB, sVSR2, IDSTT c==> Output: ULR, partial derivatives dULRdCm & radial derivative dULRdR c----------------------------------------------------------------------- c** Original Version from Nike Dattani in June 2011 for 3x3 case c** Generalized to incorporate 2x2 case, removed retardation terms and c incorporate damping ... by Kai Slaughter: July 2014 c* rj: C6{adj} & C9{adj} included in CmEFF & fixed dampF call Jan 2016 c----------------------------------------------------------------------- INTEGER NCMMax c----------------------------------------------------------------------- REAL*8 RDIST,VLIM,Cm(NCMMax),ULR,dULRdCm(NCMMax),dULRdR,R2,R3,R5, 1 R6,R8,R9,T1,T0,T2,T0P,T0P23,Dm(NCMMax),Dmp(NCMMax), 2 Dmpp(NCMMax),rhoAB,A(3,3),DR(3,3),Q(3,3),DMx(NCMMax,3,3), 3 DMtemp(3,3),DEIGMx(NCMMax,1,1),DEIGMtemp(1,1),DEIGR(1,1), 4 EIGVEC(3,1),RESID(3,1),W(3),RPOW(NCMMax),DELTAE,Modulus,Z INTEGER H,I,J,K,L,M,X,NCMM,MMLR(NCMMax),sVSR2,IDSTT,MMtemp c----------------------------------------------------------------------- DELTAE=Cm(1) R2= 1.d0/RDIST**2 R3= R2/RDIST R5= R2*R3 R6= R3*R3 R8= R6*R2 c----------------------------------------------------------------------- c....... for rhoAB.le.0.0 returns Dm(m)=1 & Dmp(m)=Dmpp(m)=0 CALL dampF(RDIST,rhoAB,NCMM,NCMMAX,MMLR,sVSR2,IDSTT,Dm,Dmp,Dmpp) c----------------------------------------------------------------------- IF(MMLR(1).GE.-1) THEN !! For the A (0) or b (-1) state c*********************************************************************** c************* Aubert Frecon 2x2 case NCMM= 7 and ... c*** Cm(1) = DELTAE c*** Cm(2) = C3Sig c*** Cm(3) = C3Pi c*** Cm(4) = C6Sig c*** Cm(5) = C6Pi c*** Cm(6) = C8Sig c*** Cm(7) = C8Pi c*********************************************************************** T1= R3*(Dm(2)*(Cm(2)-Cm(3)) + R3*Dm(4)*(Cm(4)-Cm(5)) + 1 R5*Dm(6)*(Cm(6)-Cm(7)))/3.d0 T0= DSQRT((T1 - Cm(1))**2 + 8.d0*T1**2) ULR= 0.5d0*(-Cm(1) + R3*(Dm(2)*(Cm(2)+Cm(3)) + 1 R3*Dm(4)*(Cm(4)+Cm(5)) + R5*Dm(6)*(Cm(6)+Cm(7))) + T0) c----------------------------------------------------------------------- IF(MMLR(1).EQ.0) THEN ULR= ULR + Cm(8)*R3*R6 !! add C9{adj correction ENDIF c... adjustment for the b-state IF(MMLR(1).EQ.-1) THEN ULR=ULR-T0 ULR= ULR + Cm(9)*R3*R6 !! add C9{adj correction ENDIF c... now get derivatives T0P= 0.5d0*(9.d0*T1 - Cm(1))/T0 T0P23= 0.5d0 + T0P/3.d0 c... another adjustment for the b-state IF(MMLR(1).EQ.-1) T0P23=T0P23-2.d0*T0P/3.d0 dULRdCm(1)= 0.d0 dULRdCm(2)= R3*(T0P23) dULRdCm(3)= R3*(1.d0-T0P23) dULRdCm(4)= R6*(T0P23) dULRdCm(5)= R6*(1.d0 - T0P23) dULRdCm(6)= R8*T0P23 dULRdCm(7)= R8*(1.d0-T0P23) T2 =-T0P*R3*((Dm(2)*(Cm(2)-Cm(3))+R3*(Dm(4)*2.d0*(Cm(4) 1 -Cm(5))+R2*Dm(6)*8.d0/3.d0*(Cm(6)-Cm(7))))/RDIST 2 +(Dmp(2)*(Cm(2)-Cm(3))+R3*Dmp(4)*(Cm(4)-Cm(5))+ 3 R2*R3*Dmp(6)*(Cm(6)-Cm(7)))/3.d0) dULRdR = -R3*((1.5d0*Dm(2)*(Cm(2)+Cm(3)) + R3*(Dm(4)*3.d0* 1 (Cm(4)+Cm(5))+4.d0*Dm(6)*R2*(Cm(6)+Cm(7))))/RDIST 2 + 0.5d0*(Dmp(2)*(Cm(2)+Cm(3)) + Dmp(4)*R3*(Cm(4)+ 3 Cm(5)) + Dmp(6)*R3*R2*(Cm(6)+Cm(7)))) + T2 c... and a final adjustment for the b-state IF(MMLR(1).EQ.-1) dULRdR= dULRdR- 2.d0*T2 c----------------------------------------------------------------------- ELSE c*********************************************************************** c********* Aubert Frecon 3x3 case NCMM= 10 and ... c********* Cm(1) = DELTAE c********* Cm(2) = C3Sig c********* Cm(3) = C3Pi1 c********* Cm(4) = C3Pi3 c********* Cm(5) = C6Sig c********* Cm(6) = C6Pi1 c********* Cm(7) = C6Pi3 c********* Cm(8) = C8Sig c********* Cm(9) = C8Pi1 c********* Cm(10)= C8Pi3 c*********************************************************************** c... Initialize interaction matrix to 0.d0 DO I= 1,3 DO J= 1,3 A(I,J)=0.0D0 DR(I,J)=0.d0 DO K= 1,NCMMax DMx(K,I,J)=0.d0 ENDDO ENDDO ENDDO c... Prepare interaction matrix A DO I= 2,NCMM,3 RPOW(I)= RDIST**MMLR(I) A(1,1)= A(1,1) - Dm(I)*(Cm(I)+Cm(I+1)+Cm(I+2))/(3.d0*RPOW(I)) A(1,2)= A(1,2) - Dm(I)*(Cm(I+2)+Cm(I+1)-2.d0*Cm(I))/(RPOW(I)) A(1,3)= A(1,3) - Dm(I)*(Cm(I+2)-Cm(I+1))/(RPOW(I)) A(2,2)= A(2,2) - Dm(I)*(Cm(I+2)+Cm(I+1)+4.d0*Cm(I)) 1 /(6.d0*RPOW(I)) A(3,3)= A(3,3) - Dm(I)*(Cm(I+2)+Cm(I+1))/(2.d0*RPOW(I)) ENDDO A(1,1) = A(1,1) + VLIM A(1,2) = A(1,2)/(3.d0*DSQRT(2.d0)) A(2,1) = A(1,2) A(2,2) = A(2,2) + VLIM + DELTAE A(2,3) = A(1,3)/(2.d0*DSQRT(3.d0)) A(1,3) = A(1,3)/(DSQRT(6.d0)) A(3,1) = A(1,3) A(3,2) = A(2,3) A(3,3) = A(3,3) + VLIM + DELTAE c... Prepare radial derivative of interaction matrix (? is it needed ?) DO I= 2,NCMM,3 DR(1,1)= DR(1,1) + Dm(I)*MMLR(I)*(Cm(I)+Cm(I+1)+Cm(I+2)) 1 /(3.d0*RPOW(I)*RDIST) 2 -Dmp(I)*(Cm(I)+Cm(I+1)+Cm(I+2))/(3.d0*RPOW(I)) DR(1,2)= DR(1,2) + Dm(I)*MMLR(I)*(Cm(I+2)+Cm(I+1)-2.d0* 1 Cm(I))/(RPOW(I)*RDIST) 2 -Dmp(I)*(Cm(I+2)+Cm(I+1)-2.d0*Cm(I))/(RPOW(I)) DR(1,3)= DR(1,3) + Dm(I)*MMLR(I)*(Cm(I+2)-Cm(I+1)) 1 /(RPOW(I)*RDIST) 2 -Dmp(I)*(Cm(I+2)-Cm(I+1))/(RPOW(I)) DR(2,2)= DR(2,2) + Dm(I)*MMLR(I)*(Cm(I+2)+Cm(I+1)+ 1 4.d0*Cm(I))/(6.d0*RPOW(I)*RDIST) 2 -Dmp(I)*(Cm(I+2)+Cm(I+1)+4.d0*Cm(I)) 3 /(6.d0*RPOW(I)) DR(3,3)= DR(3,3) + Dm(I)*MMLR(I)*(Cm(I+2)+Cm(I+1)) 1 /(2.d0*RPOW(I)*RDIST) 2 -Dmp(I)*(Cm(I+2)+Cm(I+1))/(2.d0*RPOW(I)) ENDDO DR(1,2) = DR(1,2)/(3.d0*DSQRT(2.d0)) DR(2,1) = DR(1,2) DR(2,3) = DR(1,3)/(2.d0*DSQRT(3.d0)) DR(1,3) = DR(1,3)/(DSQRT(6.d0)) DR(3,1) = DR(1,3) DR(3,2) = DR(2,3) c... Partial derivatives of interaction matrix A w.r.t. Cm's DO I= 2,NCMM,3 DMx(I,1,1)= -Dm(I)/(3.d0*RPOW(I)) DMx(I+1,1,1)= DMx(I,1,1) DMx(I+2,1,1)= DMx(I,1,1) DMx(I,1,2)= 2.d0*Dm(I)/(3.d0*DSQRT(2.d0)*RPOW(I)) DMx(I+1,1,2)= -DMx(I,1,2)/2.d0 DMx(I+2,1,2)= DMx(I+1,1,2) DMx(I,2,1)= DMx(I,1,2) DMx(I+1,2,1)= DMx(I+1,1,2) DMx(I+2,2,1)= DMx(I+2,1,2) DMx(I,1,3)= 0.d0 DMx(I,3,1)= 0.d0 DMx(I+1,1,3)= Dm(I)/(DSQRT(6.d0)*RPOW(I)) DMx(I+1,3,1)= DMx(I+1,1,3) DMx(I+2,1,3)= -DMx(I+1,1,3) DMx(I+2,3,1)= DMx(I+2,1,3) DMx(I,2,2)= 2.d0*Dm(I)/(3.d0*RPOW(I)) DMx(I+1,2,2)= DMx(I,2,2)/4.d0 DMx(I+2,2,2)= DMx(I+1,2,2) DMx(I,2,3)= 0.d0 DMx(I,3,2)= 0.d0 DMx(I+1,2,3)= Dm(I)/(2.d0*DSQRT(3.d0)*RPOW(I)) DMx(I+1,3,2)= DMx(I+1,2,3) DMx(I+2,2,3)= -DMx(I+1,2,3) DMx(I+2,3,2)= DMx(I+2,2,3) DMx(I,3,3)= 0.d0 DMx(I+1,3,3)= Dm(I)/(2.d0*RPOW(I)) DMx(I+2,3,3)= DMx(I+1,3,3) ENDDO c... Call subroutine to prepare and invert interaction matrix A CALL DSYEVJ3(A,Q,W) L=1 c... Now - identify the lowest eigenvalue of A and label it L DO J=2,3 IF (W(J) .LT. W(L)) THEN L=J ENDIF ENDDO c... Identifiy the highest eigenvalue of A and label it H H=1 DO J=2,3 IF(W(J).GT.W(H)) THEN H=J ENDIF ENDDO c... Identify the middle eigenvalue of A and label it M M=1 DO J=2,3 IF((J.NE.L).AND.(J.NE.H)) M= J ENDDO c... Select which eigenvalue to use based on user input IF(MMLR(1).EQ.-2) THEN X = L ELSEIF(MMLR(1).EQ.-3) THEN X = M ELSE X = H ENDIF c... determine ULR and eigenvectors ULR= -W(X) IF(MMLR(1).EQ.-2) ULR= ULR+ Cm(11)*R3*R6 !! C9adj term IF((MMLR(1).EQ.-3).OR.(MMLR(1).EQ.-4)) ULR = ULR + DELTAE IF(MMLR(1).EQ.-3) ULR= ULR+ Cm(12)*R3*R6 !! C9adj term IF(MMLR(1).EQ.-4) ULR= ULR+ Cm(13)*R3*R6 !! C9adj term DO I=1,3 EIGVEC(I,1) = Q(I,X) ENDDO cc loop over values of m to determine partial derivatives w.r.t. each Cm DO I=2,NCMM DMtemp(1:3,1:3) = DMx(I,1:3,1:3) DEIGMtemp= -MATMUL(TRANSPOSE(EIGVEC),MATMUL(DMtemp,EIGVEC)) dULRdCm(I)= DEIGMtemp(1,1) ENDDO DEIGR = -MATMUL(TRANSPOSE(EIGVEC),MATMUL(DR,EIGVEC)) dULRdR= DEIGR(1,1) !! radial derivative w.r.t. r (I think!) c------------------------------------------------------------------------ ENDIF c------------------------------------------------------------------------ RETURN CONTAINS c======================================================================= SUBROUTINE DSYEVJ3(A, Q, W) c ---------------------------------------------------------------------- c** Subroutine to setup and diagonalize the matrix A and return c eigenvalues W and eigenvector matrix Q INTEGER N, I, X, Y, R PARAMETER (N=3) REAL*8 A(3,3), Q(3,3), W(3) REAL*8 SD, SO, S, C, T, G, H, Z, THETA, THRESH c Initialize Q to the identitity matrix c --- This loop can be omitted if only the eigenvalues are desired --- DO X = 1, N Q(X,X) = 1.0D0 DO Y = 1, X-1 Q(X, Y) = 0.0D0 Q(Y, X) = 0.0D0 ENDDO ENDDO c Initialize W to diag(A) DO X = 1, N W(X) = A(X, X) ENDDO c Calculate SQR(tr(A)) SD = 0.0D0 DO X = 1, N SD = SD + ABS(W(X)) ENDDO SD = SD**2 c Main iteration loop DO 40 I = 1, 50 c Test for convergence SO = 0.0D0 DO X = 1, N DO Y = X+1, N SO = SO + ABS(A(X, Y)) ENDDO ENDDO IF(SO .EQ. 0.0D0) RETURN IF(I .LT. 4) THEN THRESH = 0.2D0 * SO / N**2 ELSE THRESH = 0.0D0 END IF c Do sweep DO 60 X = 1, N DO 61 Y = X+1, N G = 100.0D0 * ( ABS(A(X, Y)) ) IF ( I .GT. 4 .AND. ABS(W(X)) + G .EQ. ABS(W(X)) 1 .AND. ABS(W(Y)) + G .EQ. ABS(W(Y))) THEN A(X, Y) = 0.0D0 ELSE IF (ABS(A(X, Y)) .GT. THRESH) THEN c Calculate Jacobi transformation H = W(Y) - W(X) IF ( ABS(H) + G .EQ. ABS(H) ) THEN T = A(X, Y) / H ELSE THETA = 0.5D0 * H / A(X, Y) IF (THETA .LT. 0.0D0) THEN T= -1.0D0/(SQRT(1.0D0 + THETA**2)-THETA) ELSE T= 1.0D0/(SQRT(1.0D0 + THETA**2) + THETA) END IF END IF C = 1.0D0 / SQRT( 1.0D0 + T**2 ) S = T * C Z = T * A(X, Y) c Apply Jacobi transformation A(X, Y) = 0.0D0 W(X) = W(X) - Z W(Y) = W(Y) + Z DO R = 1, X-1 T = A(R, X) A(R, X) = C * T - S * A(R, Y) A(R, Y) = S * T + C * A(R, Y) ENDDO DO R = X+1, Y-1 T = A(X, R) A(X, R) = C * T - S * A(R, Y) A(R, Y) = S * T + C * A(R, Y) ENDDO DO R = Y+1, N T = A(X, R) A(X, R) = C * T - S * A(Y, R) A(Y, R) = S * T + C * A(Y, R) ENDDO c Update eigenvectors c --- This loop can be omitted if only the eigenvalues are desired --- DO R = 1, N T = Q(R, X) Q(R, X) = C * T - S * Q(R, Y) Q(R, Y) = S * T + C * Q(R, Y) ENDDO END IF 61 CONTINUE 60 CONTINUE 40 CONTINUE WRITE(6,'("DSYEVJ3: No convergence.")') END SUBROUTINE DSYEVJ3 END SUBROUTINE AFdiag c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c*********************************************************************** SUBROUTINE dampF(r,rhoAB,NCMM,NCMMAX,MMLR,sVRS2,IDSTT,DM,DMP,DMPP) c** Subroutine to generate values 'Dm' and its first `Dmp' and second c 'Dmpp' derivatives w.r.t. r of the chosen form of the damping c function, for m= 1 to MMAX. c---------------------- RJL Version of 30 January 2016 ------------------ c----------------------------------------------------------------------- c Upon Input c* r - the radial distance in Angsroms (!) c* RHOab 'universal' scaling coefficient used for systems other than H_2 c RHOab= 2*(RHOa*RHOb)/(RHOa+RHOb) where RHOa = (I_p^A/I_p^H)^0.66 c where I_p^A is the ionization potential of atom A c and I_p^H is the ionization potential of atomic hydrogen c* NCMM the number of inverse-power terms to be considered c* MMLR are the powers of the NCMM inverse-power terms c* sVRS2 defines damping s.th. Dm(r)/r^m --> r^{sVRS2/2} as r --> 0 c* IDSTT specifies damping function type: > 0 use Douketis et al. form c if IDSTT .LE. 0 use Tang-Toennies form c----------------------------------------------------------------------- c Upon Output c DM(m) - The value of the damping function for the long range term c C_MMLR(m)/r^MMLR(m) {m= 1, NCMM} c DMP(m): first derivative of the damping function DM(m) w.r.t. r c DMPP(m): second derivative of the damping function DM(m) w.r.t. r c IF(rhoAB.LE.0.0) return w. DM(m)= 1.0 & DMP(m)=DMPP(m)=0.0 for all m c----------------------------------------------------------------------- INTEGER NCMM,NCMMAX,MMLR(NCMMAX),sVRS2,IDSTT,sVRS2F,FIRST, Lsr,m, 1 MM,MMAX,MMTEMP REAL*8 r,rhoAB,bTT(-2:2),cDS(-4:4),bDS(-4:4),aTT,br,XP,YP, 1 TK, DM(NCMMAX),DMP(NCMMAX),DMPP(NCMMAX),SM(-3:25), 2 bpm(20,-4:0), cpm(20,-4:0),ZK c------------------------------------------------------------------------ c The following values for the numerical factors used in both TT and DS c were normalized to the Hydrogen data presented c by Kreek and Meath in J.Chem.Phys. 50, 2289 (1969). c The ratio has been chosen such that b= FACTOR*(I_p^X / I_p^H)^{2/3} c for the homoatomic diatomic species X_2, where I_p^A is the ionization c------------------------------------------------------------------------ DATA bTT/2.10d0,2.44d0,2.78d0,3.13d0,3.47d0/ DATA bDS/2.50d0,2.90d0,3.30d0,3.69d0,3.95d0,0.d0,4.53d0,0.d0, 1 4.99d0/ DATA cDS/0.468d0,0.446d0,0.423d0,0.405d0,0.390d0,0.d0,0.360d0, 1 0.d0,0.340d0/ DATA FIRST/ 1/ SAVE FIRST, bpm, cpm c----------------------------------------------------------------------- MMTEMP = MMLR(1) IF(MMLR(1).LE.0) MMLR(1) = 1 IF(RHOab.LE.0) THEN DO m=1,NCMMax DM(m)=1.d0 DMP(m)= 0.d0 DMPP(m)= 0.d0 ENDDO RETURN ENDIF IF(IDSTT.LE.0) THEN c=========================================== c** For Tang-Toennies type damping functions c=========================================== Lsr= sVRS2/2 IF((sVRS2.LT.-4).OR.(sVRS2.GT.4).OR.((2*LSR).NE.sVRS2)) THEN WRITE(6,600) 'TT',sVRS2 STOP ENDIF MMAX= MMLR(NCMM) + Lsr - 1 aTT= RHOab*bTT(Lsr) br= aTT*r XP= DEXP(-br) SM(-3)= 0.d0 SM(-2)= 0.d0 SM(-1)= 0.d0 SM(0)= 1.d0 TK= 1.d0 IF(br.GT.0.5d0) THEN DO m= 1,MMAX TK= TK*br/DFLOAT(m) SM(m)= SM(m-1)+ TK ENDDO DO m= 1, NCMM MM= MMLR(m) - 1 + Lsr DM(m)= 1.d0 - XP*SM(MM) DMP(m)= aTT*XP*(SM(MM) - SM(MM-1)) DMPP(m)= -aTT*aTT*XP*(SM(MM) 1 - 2.d0*SM(MM-1) + SM(MM-2)) ENDDO c----------------------------------------------------------------------- c The above section handles the calculation of the value of the damping c function for most values of r. However, at very small r that algorithm c becomes unstable due to numerical noise. To avoid this, if the c argument is very small it is re-evaluated as a finite sum ... c----------------------------------------------------------------------- ELSE MMAX= MMAX+5 DO m= 1, MMAX c... NOTE that here SM(m) is the m'th term (b*r)^m/m! [not a sum] SM(m)= SM(m-1)*br/DFLOAT(m) ENDDO DO m= 1, NCMM MM= MMLR(m) + Lsr DM(m)= XP*(SM(MM)+ SM(MM+1)+ SM(MM+2)+ SM(MM+3) 1 + SM(MM+4)) DMP(m)= aTT*XP*SM(m-1) DMPP(m)= aTT*aTT*XP*(SM(m-2)-SM(m-1)) ENDDO ENDIF ENDIF c IF(IDSTT.GT.0) THEN c======================================================================= c** For Douketis-Scoles-Marchetti-Zen-Thakkar type damping function ... c======================================================================= IF((sVRS2.LT.-4).OR.(sVRS2.GT.4).OR.(sVRS2.EQ.1).OR. 1 (sVRS2.EQ.3)) THEN WRITE(6,600) 'DS',sVRS2 STOP ENDIF IF(FIRST.EQ.1) THEN DO m= 1, 20 DO sVRS2F= -4,0 bpm(m,sVRS2F)= bDS(sVRS2F)/DFLOAT(m) cpm(m,sVRS2F)= cDS(sVRS2F)/DSQRT(DFLOAT(m)) ENDDO ENDDO FIRST= 0 ENDIF br= rhoAB*r DO m= 1, NCMM MM= MMLR(m) XP= DEXP(-(bpm(MM,sVRS2) + cpm(MM,sVRS2)*br)*br) YP= 1.d0 - XP ZK= MM + 0.5d0*sVRS2 DM(m)= YP**ZK TK= (bpm(MM,sVRS2) + 2.d0*cpm(MM,sVRS2)*br)*rhoAB DMP(m) = ZK*XP*TK*DM(m)/YP c ... calculate second derivative [for DELR case] {check this!} DMPP(m)= (ZK-1.d0)*DMP(m)*(XP*TK)/YP 1 - DMP(m)*TK + DMP(m)*2.d0*cpm(MM,sVRS2)*rhoAB**2/TK ENDDO ENDIF MMLR(1) = MMTEMP RETURN 600 FORMAT(/,' *** ERROR *** For ',A2,'-damping functions not yet de 1fined for sVRS2=',i3) END c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12