C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C STRFIT - GENERAL STRUCTURE FITTING PROGRAM, FITTING EFFECTIVE STRUCTURAL C PARAMETERS DEFINED IN THE 'CART' SCHEME TO ANY COMBINATION C OF ROTATIONAL CONSTANTS C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C This program was originally written in 1990 for determination of C effective r.0 structures of hydrogen bonded complexes, but currently it C can also be used to fit also r.z, rm(2L) and other geometries C C Citation and additonal information on this program: C Z. Kisiel, J. Mol. Spectrosc. 218, 58-67 (2003). C C This STRFIT was written independently of, but seems to embody similar C features to Schwendeman's STRFIT. Once I found this out it was too late C to change the name. C C See: R.H.Schwendeman, "Structural parameters from rotational spectra", C in "Critical evaluation of chemical and physical structural C information", D.R.Lide and M.A.Paul, Eds., National Academy of C Science, Washingtom, D.C. (1974). C C for a similar least-squares fitting approach see also: C C P.Nosberger,A.Bauder, and Hs.H.Gunthard, "A versatile method C for molecular structure determinations from ground state C rotational constants", Chem.Phys. 1,418-425(1973). C C The least-squares engine for fitting internal parameters to C moments of inertia can now be used to evaluate also: C C 1/ r.z geometry or any other geometry fitted to moments of inertia C modified by user supplied corrective contributions such as C r_se geometry (=semi empirical based on B0's + ab initio alphas) C 2/ r.m(1r) and rm(1rL) geometry C 3/ r.m(2r) and rm(2rL) geometry C C For r.m(1), r.m(2) and their subvariants see: C C J.K.G.Watson, A.Roytburg, W.Ulrich, J.Mol.Spectrosc. C 196,102-119(1999) C C Declaration scheme of internal structural coordinates C (colloquially called CART since it was adapted from a Wilson C group program of that name): C C H.R.Thompson, J.Chem.Phys. 47,3407-3410(1967) C C C Main features of STRFIT: C C - fit of up to MAXPAR structural parameters to up to MAXCON measured C rotational constants in up to MAXISO isotopic species C - each structural parameter can have up to MAXREP repetitions in the C molecule C - up to MAXCHG structural parameters can be different between the parent C and the substituted species C - non-linear Levenberg-Marquardt least squares fitting algorithm is used C (Numerical Recipes Chapt.14) C - the fit is to (equally weighted) moments of inertia or planar moments C and not to rotational constants, although when rotational constants C are to be fitted they are declared conventionally C - the spectroscopic constants that can be fitted to are A,B,C, B+C, C and A+B, and the three planar moments Pa,Pb,Pc (in any combination) C - the data set for STRFIT is compatible with PMIFST so the latter can C be used to check the 'CART' definitions C - the results of STRIFT are written to the main output file and also to C file STRFIT.PMI for direct inspection with PMIFST C C C ver. 15.II.2010 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C C 1990: Creation C 24.06.95: PARAMETER dimensioning C 10.12.99: Fitting of B+C and A+B C 30.08.00: Echo of fitted structure to STRFIT.PMI C 8.01.01: Various incremental modifications C 8.04.02: Fitting also of planar moments C 19.04.02: Elimination of bug in writing STRFIT.PMI C 4.06.02: rm(1) and rm(2) model + fixed parameters + overhaul C 29.10.02: increase in MAXPAR and associated corrections C 30.01.03: fitting equal c or d constants C 8.08.05: correction of 'final principals' bug spotted by Stew Novick C 14.11.08: correction to mass used for d terms, spotted by Jean Demaison C 25.05.09: elimination of erroneous temporary code that crept in about 2005 C in Laurie parameter + improvement of symmetric top performance C 15.02.10: allowance also for dB=B*-Bexp input and better input debugging C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C DATA FILE STRUCTURE - specimen data file is reproduced below, with C lines numbered to facilitate description: C C first column of the data file C | 32nd column of data file (beginning of C | | numerical input for all lines C | | below the CART deck) C lno. | | C 1. buteneyne...HF, HF at the TRIPLE bond, in the bey plane C 2. 11 0 C 3. 1 0 0 0 .000000 .000000 .000000 12.000000 C 4. 2 1 0 0 1.208600 .000000 .000000 12.000000 C 5. 3 2 1 0 1.431000 177.900000 .000000 12.000000 C 6. 4 3 2 1 1.341600 123.100000 180.000000 12.000000 C 7. 5 1 2 4 1.062000 182.300000 180.000000 1.007825 C 8. 6 4 3 2 1.087000 121.600000 .000000 1.007825 C 9. 7 4 3 2 1.087000 118.700000 180.000000 1.007825 C 10. 8 3 4 7 1.087000 121.700000 0.000000 1.007825 C 11. 9 2 3 8 0.604300 177.900000 .000000 .000001 C 12. 10 9 2 3 3.080000 90.000000 180.000000 18.998404 C 13. 11 10 9 2 0.925700 360.000000 .000000 1.007825 C 14. TOTAL NUMBER OF PARAMETERS: 2 C 15. atom no., parameter no. 10 1 0 C 16. atom no., parameter no. 9 1 0 C 17. NO OF CONSTANTS TO FIT TO: 6 C 18. constant,species,value 1 1 5182.82917 C 19. constant,species,value 2 1 3312.90406 C 20. constant,species,value 3 1 2021.475817 C 21. constant,species,value 1 2 5182.55955 C 22. constant,species,value 2 2 3271.83826 C 23. constant,species,value 3 2 2005.91671 C 24. NO OF CHANGES FROM PARENT SP.: 2 C 25. atom no.,parameter no.,value 11 4 2.0141022 C 26. atom no.,parameter no.,value 11 1 -0.0023 C C where: C C line.1 - comment C line.2 - the number of atoms in the molecule (read in with I5) C optionally a second number (also in I5 field) can be included, C and if equal to 1 specifies output of various debugging ] C information C line.3-13 - CART definition of the molecule - see a separate set of comments C further below. Defined centres can be either real or dummy atoms C and the use of dummy atoms is encouraged since it often C simplifies things considerably. In the example above a dummy atom C (carrying an infinitesimally small but non-zero mass) defines C the centre of the triple bond to which the HF molecule C is notionally bonded. The program will actually use the mass in C the calculation so specify as low as possible (but non-zero mass) C and mass of 0.000001u is recommended. C Dummy atoms should not be counted in the evaluation of the rm^2 C d-parameters and are identified as such if their mass is less C than 0.00001u. C line.14+ - these lines are read in with (a31,2i2,..), the contents of the C alphanumeric field is for descriptive purposes only. C The lines contain: definitions of parameters to C be fitted, observed values of rotational constants and, if C required, definitions of changes relative to the parent C species if constants are for more than one species. C C line.14: The total number of parameters, being a total of parameters to C be fitted and of parameters to be held fixed at supplied values. C The parameters can be either internal structural coordinates, or C ancillary parameters such as c,d parameters of the rm(1) and C rm(2) models, Laurie delta_H parameters, or empirical deuteration C bond length change C C line.15: DEFINITION OF AN INTERNAL COORDINATE C ------------------------------------ C - the first number is the atom number, C - the second number is the CART parameter for that atom C (1=bond length, 2=bond angle, 3=dihedral angle) C - the third number is an optional repeat count and defines whether C this parameter does not have additional occurrences: C a/ If the parameter is unique then set the repeat count to zero C or ignore it altogether. C b/ Alternatively a non-zero repeat count makes it is possible C to fit 'degenerate' parameters, occurring for more than one C atom (eg. CX bond in CX3 pyramidal segment): in such case C the third parameter in this line should be set to the number C of additional occurrences of the parameter (e.g. 2 for bond C lengths in a methyl group), and an additional line C (in a31,20i2 format) should be inserted C immediately following, with the numbers of atoms also carrying C this parameter. C Up to ten such additional occurrences are allowed. C NOTE: symmetric +- dihedral angle values can be fitted by C specifying the values which have opposite sign to that of C the parent value by using a negative number for the respective C atom. Regrettably, owing to current formatting restrictions, C the numbers of such of such atoms have to be smaller C than 10 (single digit) C - the optional fourth number is an initial value for this parameter, C which will override the value in the preceding CART deck C C C line.15a: DEFINITION OF AN ANCILLARY PARAMETER OF FIT C --------------------------------------------------- C Example (line read in with format 31x,3i2,f14.0): C C rm(1) parameter................-4 1 0 0.001 C C - the first number selects the type of the ancillary parameter C = -4 Watson's rm(1) c parameter, in which case the second number C = 1,2,3,4,5,6 for c_a,c_b,c_c,c_ab,c_ac,c_bc C = 7 for fitting of c_a = c_b C = 8 for fitting of c_a = c_c C = 9 for fitting of c_b = c_c C =10 for fitting c_a = c_b = c_c C = -5 Watson's rm(2) d parameter, in which case the second number C = 1,2,3 for d_a,d_b,d_c C = 7 for fitting of d_a = d_b C = 8 for fitting of d_a = d_c C = 9 for fitting of d_b = d_c C =10 for fitting d_a = d_b = d_c C = -6 Laurie correction term for an XH bond, in which case C the second number specifies the atom for which this bond C is defined C = -7 Laurie term defined in terms of a projection angle C = -8 isotopic bond length change relative to the value in C the parent - this is only for use with Deuterium C substitution. The second number is then the number of the C atom for which this bond has been defined C C - the third number: C - for -4,-5 = 0 C - for -6,-7,-8 = specifies the number of additional C internals subject to the same correction and an additional C line (in a31,20i2 format) should be inserted immediately C following, with numbers of atoms for which these internals C are defined C C - the fourth (floating point) number is the initial value C for the parameter C C NOTE C ---- C Any of the declared internal or ancillary parameters can be fixed C instead of fitted by writing FIX in columns 29-31 of the C declaration line - this is the only case for which the text in C in columns 1-31 is relevant. Otherwise this is ignored and C can be used for explanatory comments C C C line.18: DEFINITION OF THE SPECTROSCOPIC CONSTANT TO BE FITTED TO: C --------------------------------------------------------- C line read in with format (31x,2i2,3f14.0): C ++--++++++++++++++--------------++++++++++++++ C Examples: C C Constant,species,value....... 1 1 312142.390 NH2CN C Constant,species,value....... 1 1 312142.390 1.0 NH2CN C Constant,species,value....... 1 1 312142.390 30000. NH2CN C C - the first number defines the type of constant C = 1,2,3 for A,B,C C = 4,5 for B+C and A+B C = 6,7,8 for P.a,P.b,P.c, C - the second number defines the isotopic species, C - the third number is the value of the constant, in MHz for C rotational constants and uA**2 for planar moments. C The constants should be declared in C increasing order of isotopic species number, all constants for C a given species are not required. c c It is possible to specify optional fourth OR fifth numbers c (only one allowed), such that: C C - the fourth number is a vibration-rotation correction C dI=Icorr-Iobs (in uA**2) C to the experimental moment of inertia defined by the preceding C third number. C - the fifth number is a vibration-rotation correction C dB=Bcorr-Bobs (in MHz) C to the experimental rotational constant defined by the preceding C third number. C C The structure will then be fitted to corrected quatnities C Icorr=Iobs+dI, Iobs=conv/Bobs or to C Icorr=conv/(Bobs+dB) C where Iobs is the moment of inertia derived from the preceding rotational C constant, or to P=Pobs+dI if planar moment has been specified C C --NOTE: the dI option is built in for use with various models of C subtracting vibration-rotation contributions such as C a/ r* (r.z) structures by using dI's from harmonic force field C b/ semi-experimental structures by using dB's from anaharmonic C force field C --NOTE: this line is read in fixed format, and the value of C constant and optional dI are read in with 2F14.x C --NOTE: when B+C and A+B are declared, the fit is actually to C 4*505379.01/(B+C) = 4*I.bI.c/(I.b+I.c) =ca. I.b+I.c or to C 4*505379.01/(A+B) = 4*I.aI.b/(I.a+I.b) =ca. I.a+I.b C C line.24-26 - DEFINITION OF ISOTOPIC CHANGES C ------------------------------ c These lines, and more, if necessary only follow if rotational C constants are for more than one isotopic species. These sets C of lines should be placed in the order in which the isotopic C species are numbered in lines declaring the rotational C constants. C line.25: definition of changes relative to the parent species C - the first number refers to the atom C - the second number is the parameter (1,2,3,4 for C bond length, bond angle, dihedral angle and atomic mass resp.), C - the third number is the value of this parameter in C the substituted species. When this parameter is not a mass, C the number is the change in the parameter relative to the C parent species C C GENERAL NOTES: C C - As noted in several places above the input is in the C old-fashioned fixed format, so it is necessary to keep input C numbers within the specified columns. C C - Columns to the right of the input fields can be used for comments. C In most cases this applies to columns 65 and higher, but for C input of spectroscopic constants such columns are 78 and higher. C C - The input file is not read past the last set of isotopic changes, C so the rest of the file can be used for any purpose - usually C storage of alternative data and of results. C C - STRFIT does not write to the input file so it will not corrupt C it - but can overwrite it if its name is accidentally given C as the output file, so take care C C - Additional blocks of comments are allowed before all lines C with a repeat count: ie. those that specify: C the number of atoms, C the number of parameters, C the number of spectroscopic constants C the number of isotopic changes. C These lines have to have a '!' character in the first column, C and there is no limit on their number C C WARNING: There is a known problem with symmetric top molecules C which break symmetry on isotopic substitution. This is C associated with the way the program implements parallel axes for C the rm(1r) scheme and possibility of perpendicular axis switching. C Steps taken to stabilise this performance are marked by the string C May2009 C C----------------------------------------------------------------------------- C C O M P I L A T I O N: C----------------------------------------------------------------------------- C C If you want to compile this program then please note that it is crucially C dependent on the use of 'static memory allocation'. Please read the following C notes copied from the PROSPE webpage (and also use the files from the webpage C for testing): C C Most contemporary compilers no longer use default static memory allocation C that preserves values of variables previously set in some subroutine. Many C programs assume this to be the case so that you need to use an appropriate C keyword to enforce static allocation. For example on the f77 compiler on C SGI machines this is the -static option, whereas with Intel Visual Fortran C you have to use option -Qsave. Note that with some compilers optimization C options used not to be safe. With good compilers this is no longer the C case, but if problems crop up it is still good to check by disabling C optimization. C C----------------------------------------------------------------------------- C CART = Thompson J.Chem.Phys. 47,3407(1967) definitions: C----------------------------------------------------------------------------- C C 1. Atom N in the molecule is defined in terms of three other neighbouring C atoms: NA, NB, and NC, and a bond length R(N-NA), a bond angle C A(N-NA-NB) and a dihedral angle D(N-NA-NB-NC) C C 1a. Atoms are declared in sequence of increasing N (starting from N=1) C and the numbers of the defining atoms NA, NB, NC all have to be smaller C than N (and all different) C C 1b. For N=1 set NA,NB,NC, R,A,D to zero C For N=2 set NB,NC, A,D to zero C For N=3 set NC, D to zero C C 2. For a right handed coordinate system dihedral angles are defined C by viewing the configuration C i/ from the NC direction along the NB-NA axis if NC is bonded to NB C ii/ from the NC direction along the NA-NB axis if NC is bonded to NA C The angle is defined positive for a clockwise rotation C of atom NC into the plane N-NA-NB. C C 3. In the calculation of initial cartesian coordinates the atom N=1 C defines the origin, the line N=1 to N=2 defines the X axis and C atoms N=1, N=2, N=3 define the XY plane. C C____________________________________________________________________________ C C IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL LMUSER C parameter (maxpar=50,maxiso=30,maxcon=3*maxiso,maxrep=20, * maxchg=10,maxato=300) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /INPUT/CARTRF(MAXATO,4),NPARF(maxpar,3), * NREP(maxpar,maxrep),VPARF(maxpar),NCOBS(MAXCON,2), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2), * VCHANG(maxiso,maxchg),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,COMENT COMMON /DERIV/DER(MAXCON,maxpar) COMMON /coords/CORDP(MAXATO,3) COMMON /COV/COVTMP(maxpar,maxpar) COMMON /EPSIL/ristr(maxcon,3),ritot(maxcon,3) REAL*8 X(MAXCON),Y(MAXCON),SIG(MAXCON),A(maxpar),ALAMDA,CHISQ, * COVAR(maxpar,maxpar),ALPHA(maxpar,maxpar),DYDA(maxpar), * ERRPAR(maxpar),AOLD(maxpar) INTEGER LISTA(maxpar) CHARACTER NAMPAR(maxpar)*20,coment*78,nameco(8)*3, * cnames(maxpar)*20 DATA nameco/' a ',' b ',' c ','b+c','a+b','P.a','P.b','P.c'/ c NCA=maxpar C C...Input of data C CALL DATAIN <------- C C C...Nonlinear L-M least-squares fit C DO 12 I=1,NOBSC X(I)=I Y(I)=VCOBS(I,1) SIG(I)=1.0d0 12 CONTINUE c NPT=NOBSC kk=0 DO 13 I=1,npar if(ifit(i).eq.0)goto 13 kk=kk+1 A(kK)=VPARF(I) AOLD(kK)=A(kK) LISTA(kk)=kK 13 CONTINUE MA=kK MFIT=kK c c... initialize c NDEGF=NPT-MA IF(NDEGF.EQ.0)NDEGF=1 ALAMDA=-1.d0 niter=-1 CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA, <------- * NCA,CHISQ,LMUSER,ALAMDA) c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c Main fitting loop c WRITE(*,'(1X/)') 1 WRITE(3,6006)niter,CHISQ,ALAMDA 6006 FORMAT(/1x,39(2H- )/' ITERATION',I5,': Chi-squared =',F16.10, 1 5X,'ALAMDA =',E9.2//' Constants (value,error,change):') SWDEV=SQRT(CHISQ/NDEGF) C i=0 DO 52 ii=1,npar if(ifit(ii).eq.0)goto 52 i=i+1 ERRPAR(I)=SQRT(COVTMP(I,I))*SWDEV WRITE(3,53)NAMPAR(II),A(I),ERRPAR(I),A(I)-AOLD(I) 53 FORMAT(1X,A,F12.5,' +-',F10.5,F16.6) AOLD(I)=A(I) 52 CONTINUE kk=0 do 55 i=1,npar if(ifit(i).eq.0)goto 55 kk=kk+1 vparf(i)=a(kk) 55 continue C WRITE(*,50)niter,CHISQ,ALAMDA 50 FORMAT(1X,'ITERATION',I5,': Chi-squared =',F11.5, 1 5X,'ALAMDA =',E9.2) OCHISQ=CHISQ CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA, <------- * NCA,CHISQ,LMUSER,ALAMDA) IF (ABS(OCHISQ-CHISQ).GT.1.D-10)GOTO 1 C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C Write final results: C 1/ to UNIT=3 main output file C 2/ to UNIT=4 echo of fitted structure to file STRFIT.PMI C (with masses of the parent isotopomer) C WRITE(*,'(1X/'' ===== FINAL RESULTS =====''/)') ALAMDA=0.0d0 CALL MRQMIN(X,Y,SIG,NPT,A,MA,LISTA,MFIT,COVAR,ALPHA, <------- * NCA,CHISQ,LMUSER,ALAMDA) C OPEN(4,FILE='strfit.pmi',STATUS='UNKNOWN',ERR=520) WRITE(4,'(A)')COMENT WRITE(4,'(I5)')NATMS C C...make sure CARTDF contains final parameters of fit and not values leftover C from last evaluation of derivatives C do 57 i=1,natms do 57 j=1,3 cartdf(i,j)=cartrf(i,j) 57 continue C kk=0 do 56 i=1,npar if(ifit(i).eq.1)then kk=kk+1 vparf(i)=a(kk) endif 56 continue c DO 120 NN=1,NPAR if(nparf(nn,1).gt.0)then CARTDF(NPARF(NN,1),NPARF(NN,2))=VPARF(NN) IF(NPARF(NN,3).NE.0)THEN DO 122 NNN=1,NPARF(NN,3) CARTDF(iabs(NREP(NN,NNN)),NPARF(NN,2))=VPARF(NN) if(nrep(nn,nnn).lt.0) * CARTDF(iabs(NREP(NN,NNN)),NPARF(NN,2))=-VPARF(NN) 122 CONTINUE ENDIF endif 120 CONTINUE C C...write current CART data to STRFIT.PMI C DO 3205 N=1,NATMS WRITE(4,3201)(NCART(N,NN),NN=1,4),(CARTDF(N,NN),NN=1,3), * CARTRF(N,4) 3205 CONTINUE 3201 FORMAT(4I4,3F13.6,F16.7) write(4,'('' |''/'' |''/ * '' |__ Generated by STRFIT from a run resulting in:'' * /)') C C...write final results: fitted constants and errors C WRITE(3,'(1X,78(1H_)//'' FINAL RESULTS OF LEAST SQUARES FIT:'' 1 /)') SWDEV=SQRT(CHISQ/NDEGF) SWDSQ=SWDEV*SWDEV n=0 DO 6004 k=1,NPAR if(ifit(k).eq.1)then n=n+1 ERRPAR(N)=SQRT(COVAR(N,N))*SWDEV WRITE(*,6032)NAMPAR(k),A(N),ERRPAR(N) if(nparf(k,3).eq.0)then WRITE(3,6032)NAMPAR(k),A(N),ERRPAR(N) WRITE(4,6032)NAMPAR(k),A(N),ERRPAR(N) else WRITE(3,6033)NAMPAR(k),A(N),ERRPAR(N), * (nrep(k,nn),nn=1,nparf(k,3)) WRITE(4,6033)NAMPAR(k),A(N),ERRPAR(N), * (nrep(k,nn),nn=1,nparf(k,3)) endif else if(nparf(k,3).eq.0)then WRITE(3,6132)NAMPAR(k),vparf(k) WRITE(4,6132)NAMPAR(k),vparf(k) else WRITE(3,6133)NAMPAR(k),vparf(k), * (nrep(k,nn),nn=1,nparf(k,3)) WRITE(4,6133)NAMPAR(k),vparf(k), * (nrep(k,nn),nn=1,nparf(k,3)) endif endif 6004 CONTINUE 6032 FORMAT(1X,A,F13.6,' +-',F10.6) 6033 FORMAT(1X,A,F13.6,' +-',F10.6,' and at N=',20I3) 6132 FORMAT(1X,A,' [',F11.6,' ] FIXED ') 6133 FORMAT(1X,A,' [',F11.6,' ] FIXED ',' and at N=',20I3) C WRITE(3,6005)CHISQ,SWDEV WRITE(4,6005)CHISQ,SWDEV WRITE(*,6005)CHISQ,SWDEV WRITE(*,'(1X/)') 6005 FORMAT(/' Chi-squared =',F20.10/ 1 ' Deviation of fit =',F16.6) c c...Obs. and Obs.-Calc. quantities c WRITE(3,23) WRITE(4,23) 23 FORMAT(/'Ni Axis',5X,'Iobs',7X,'Icalc',5X,'Io-c',13X, 1 'Bobs',7X,'Bcalc',5X,'Bo-c'/) DO 20 I=1,NPT WRITE(3,22)ncobs(i,2),nameco(ncobs(i,1)), 1 VCOBS(I,1),VCOBS(I,2),VCOBS(I,1)-VCOBS(I,2), 1 VCOBS(I,3),VCOBS(I,4),VCOBS(I,3)-VCOBS(I,4) WRITE(4,22)ncobs(i,2),nameco(ncobs(i,1)), 1 VCOBS(I,1),VCOBS(I,2),VCOBS(I,1)-VCOBS(I,2), 1 VCOBS(I,3),VCOBS(I,4),VCOBS(I,3)-VCOBS(I,4) 20 CONTINUE 22 FORMAT(i2,1x,a,2F12.5,F10.5,5X,2F12.4,F10.4) C C...correlation coefficients (code lifted from PISORT<-ASFIT) C WRITE(3,32) WRITE(4,32) 32 FORMAT(1X/' Correlation coefficients:') n=0 DO 30 k=1,NPAR if(ifit(k).eq.0)goto 30 n=n+1 DO 31 NN=1,N COVAR(N,NN)=SWDSQ*COVAR(N,NN)/(ERRPAR(N)*ERRPAR(NN)) 31 CONTINUE cnames(n)=nampar(k) 30 CONTINUE ninfit=n c M=1 305 MM=M+7 IF(MM.GT.ninfit)MM=ninfit WRITE(3,301) WRITE(3,230)(L,L=M,MM) WRITE(4,301) WRITE(4,230)(L,L=M,MM) 230 FORMAT(21x,8(i7)) WRITE(3,301) WRITE(4,301) DO 302 N=M,ninfit LLLL=MM IF(LLLL.GT.N)LLLL=N WRITE(3,303)n,cnames(N)(1:19),(covar(N,L),L=M,LLLL) WRITE(4,303)n,cnames(N)(1:19),(covar(N,L),L=M,LLLL) 302 CONTINUE 303 FORMAT(i2,':',A,8(F7.3)) 301 format(1x) IF(MM.GE.ninfit)THEN GOTO 3031 ELSE M=M+8 GOTO 305 ENDIF c 3031 WRITE(3,'(1X,78(1H_))') WRITE(4,'(1X,78(1H_))') C C...Final Cartesians to main output C do 300 i=1,natms cartdf(i,4)=cartrf(i,4) 300 continue call momin(1,0) <------- WRITE(3,3111) 3111 FORMAT(1X/' Final principal coordinates of parent:'// * 1X,'ATOM NO.',7X,'A',14X,'B',14X,'C',13X,'MASS'/) DO 3112 I=1,NATMS WRITE(3,3114)I,CORDp(I,1),CORDp(I,2),CORDp(I,3),cartrf(I,4) 3112 CONTINUE 3114 FORMAT(2X,I2,2X,3F15.6,F17.7) WRITE(3,'(1X,78(1H_))') C C...Final inertial contributions from c,d,Laurie and shift parameters C neps=0 do 6200 k=1,npar if(nparf(k,1).lt.0)neps=neps+1 6200 continue C if(neps.gt.0)then call setup1(-1,0,A,MA) <------- WRITE(3,623) 623 FORMAT(1x/' Terms in I.fitted = I.rigid + eps'/ * /' Ni',4X,'(I_a)rig',4X,'(I_b)rig',4X,'(I_c)rig',9X, * 'eps_a',7X,'eps_b^**',4X,'eps_c'/) do 6202 k=1,nisot write(3,622)k,(ristr(k,nn),nn=1,3), * (ritot(k,nn)-ristr(k,nn),nn=1,3) 6202 continue 622 FORMAT(i3,3F12.5,4X,3F12.5) WRITE(3,624) 624 FORMAT(1x/'** Only this value is the total EPSILON for the ', * 'degenerate constant in linear'/' and symmetric tops, ', * 'the other value does not contain the rm() c and d'/ * ' contributions.') WRITE(3,'(1X,78(1H_))') endif C close(3) c c...write isotopomer specifications to STRFIT.PMI c IF(NISOT.gt.1)then DO 6007 N=2,NISOT WRITE(4,6042)N DO 6007 NN=1,NUMCHG(N) WRITE(4,6043)' atom no.,parameter no.,value ', 1 (NCHANG(N,NN,I),I=1,2),VCHANG(N,NN) 6007 CONTINUE 6042 FORMAT(1X/' ISOTOPIC SPECIES',I3,', changes from ', 1 'parent species:') 6043 FORMAT(A,2I3,F17.7) WRITE(4,'(1X,78(1H_))') endif close(4) c C 520 STOP END C C_____________________________________________________________________________ C SUBROUTINE DATAIN C C............................................................................ C C INPUT OF DATA, where: C C NATMS - number of atoms in the molecule C NDATMS - number of dummy atoms in the molecule, assumed to have mass C <0.00001u (needed in conection with rm^2 d coefficients) C NCART - numbers N,NA,NB,NC in the CART definition scheme for NATMS atoms C CARTDF - bond lengths N.NA, angles N.NA.NB, dihedral angles N.NA.NB.NC and C atomic masses for NATMS atoms C CARTRF - this holds the reference (parent) CART parameters whereas CARTDF C contains parameters submitted for calculation C NPAR - number of declared constants (some of which may be fixed) C NPARF - definition of the parameter to be fitted or fixed. If this is an C internal structural coordinate then: C (i,1) - the number of atom C (i,2) - the required CARTDF parameter for that atom C (i,3) - the number of additional occurrences of this parameter in C the CART definitions C If the parameter is not an explicit internal coordinate then C (i,1) - the type of parameter C =-4 rm(1) c-parametr C =-5 rm(2) d-parameter C =-6 Laurie hydrogen atom type correction C =-7 Laurie type correction defined by projection angle C =-8 straightforward bond shrinkage C =-9 Laurie type correction for bond angle (trial only - not C tested) C (i,2) - = 1,2,3,4,5,6 for c_a,c_b,c_c,c_ab,c_ac,c_bc C 7,8,9 for c_a=c_b,c_a=c_c,c_b=c_c C 10 for c_a = c_b = c_c C = 1,2,3 for d_a,d_b,d_c C 7,8,9 for d_a=d_b,d_a=d_c,d_b=d_c C = the number of atom whose bond has bond shrinkage C (i,3) - the number of isotopic species implicated in bond shrinkage C C NREP - additional occurrences of a given parameter where: C (i,j) - is the atom number which also carries parameter i C VPARF - values of the declared parameters, if initial value for an internal C is not specified on input, it is transferred from CARTDF C IFIT - declaration of whether the parameter is fixed (=0) or to be C fitted (=1) C NOBSC - the number of measured constants available for the fit C NCOBS - definitions of the observed constants where: C (i,1) - the type of constant: = 1,2,3,4,5,6,7,8 for A,B,C,B+C,A+B, C Pa,Pb,Pc resp. C (i,2) - the number of the isotopic species, = 1 for the parent as C read into CARTDF C VCOBS - values of the observed and calculated rotational constants and C principal moments of inertia where: C (i,1) - obs. moment of inertia C (i,2) - calc. moment of inertia C (i,3) - obs. rotational constant C (i,4) - calc. rotational constant C NISOT - the number of isotopic species for which these constants have been C measured C NUMCHG - the number of parameters different for each additional species C relative to the parent one C NCHANG - coordinates of parameters changed relative to the parent species: C (i,j,k) - i refers to the number of the species C - j refers to the number of parameter being changed for that C species C - k=1 refers to number of atom carrying that parameter C - k=2 refers to the CARTDF parameter number (1,..,4) C VCHANG - values of parameters different from the parent species: C (i,j) - i is the species number C - j is the changed parameter number C............................................................................ C C NAMPAR - alphanumeric description for parameters of fit, set up in DATAIN C and used later C............................................................................ C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=50,maxiso=30,maxcon=3*maxiso,maxrep=20, * maxchg=10,maxato=300,conv=505379.01D0) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /INPUT/CARTRF(MAXATO,4),NPARF(maxpar,3), * NREP(maxpar,maxrep),VPARF(maxpar),NCOBS(MAXCON,2), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2), * VCHANG(maxiso,maxchg),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,COMENT CHARACTER CNAME(8)*3,PARNM(6)*2,dastr*26 CHARACTER*20 NAMPAR(maxpar) CHARACTER COMENT*78,line*78 CHARACTER*31 FILNAM,TEXT DATA CNAME/' A ',' B ',' C ','B+C','A+B','P.a','P.b','P.c'/ DATA PARNM/'a ','b ','c ','ab','ac','bc'/ C C...comment and CART definition scheme for the molecule - this defines the C parent species C WRITE(*,'(1x//)') WRITE(*,6002) 6002 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | STRFIT - General structure fitting', * ' program using CART definitions',T79,'|'/ * ' |',76(1H_),'|'/' version 15.II.2010',T64,'Zbigniew KISIEL'/) C 6001 WRITE(*, 3040)'Name of data file: ' 3040 FORMAT(1X/1X,A,$) READ(*,'(A)',ERR=6001)FILNAM OPEN(2,FILE=FILNAM,ERR=6001,status='old') nlin=0 C 6010 WRITE(*,3040)'Name of output file: ' READ(*,'(A)',ERR=6010)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=6010) C READ(2,'(A)',ERR=6009)LINE nlin=nlin+1 READ(LINE,'(A)',ERR=6008)COMENT write(3,6002) WRITE(3,6703)COMENT 6703 FORMAT(1x,78(1H_)//1X,A/1X,78(1H_)/) C 6702 READ(2,'(a)',err=6009)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE goto 6702 endif write(*,'(1x)') C READ(LINE,3003,ERR=6008)NATMS,idump 3003 FORMAT(2I5) DO 3206 N=1,NATMS READ(2,'(A)',ERR=6009)LINE nlin=nlin+1 READ(LINE,*,ERR=6008)(NCART(N,NN),NN=1,4),(CARTRF(N,NN),NN=1,4) 3206 CONTINUE ndatms=0 DO 6020 N=1,NATMS if(abs(cartrf(n,4)).lt.0.00001d0)ndatms=ndatms+1 DO 6020 NN=1,4 CARTDF(N,NN)=CARTRF(N,NN) 6020 CONTINUE C if(ndatms.eq.0)then dastr=' ' else write(dastr,3028)ndatms 3028 format('(including',i3,' dummy atoms)') if(ndatms.eq.1)dastr(25:26)=') ' endif c WRITE(3,3021)NATMS,dastr 3021 FORMAT(1x/ * 1X,'NUMBER OF ATOMS = ',I2,3x,a// * 1X,' NO NA NB NC NO.NA NO.NA.NB NO.NA.NB.NC', * ' MASS'/) DO 3022 N=1,NATMS WRITE(3,3023)(NCART(N,NN),NN=1,4),(CARTDF(N,NN),NN=1,4) 3022 CONTINUE 3023 FORMAT(1X,4I4,3F13.6,F15.7) C C - - - - - - - - - - - - - C C...Read parameters to fit/fix C 6700 READ(2,'(a)',ERR=6009)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE goto 6700 endif C NFIX=0 READ(line,6003,ERR=6008)TEXT,NPAR 6003 FORMAT(A,2I2,3F14.4) WRITE(3,6011)' TOTAL NUMBER OF PARAMETERS: ',NPAR 6011 FORMAT(1X/1X,A,I5) WRITE(3,'(1x/1x,a/7x,24(1h-))')' Parameters to be fitted: ' DO 6004 N=1,NPAR READ(2,'(A)',ERR=6009)LINE nlin=nlin+1 READ(LINE,5003,ERR=6008)TEXT,(NPARF(N,I),I=1,3),value 5003 FORMAT(A,3I2,f14.0) if(text(29:31).eq.'FIX')then ifit(n)=0 NFIX=NFIX+1 else ifit(n)=1 endif vparf(n)=value C C...rm() c-parameter C if(nparf(n,1).eq.-4)then I=nparf(n,2) if(i.le.6)then WRITE(NAMPAR(N),5001)'c_'//parnm(i)//' =' else if(i.eq.7) WRITE(NAMPAR(N),'(8x,a)')'c_a = c_b =' if(i.eq.8) WRITE(NAMPAR(N),'(8x,a)')'c_a = c_c =' if(i.eq.9) WRITE(NAMPAR(N),'(8x,a)')'c_b = c_c =' if(i.eq.10)WRITE(NAMPAR(N),'(2x,a)')'c_a = c_b = c_c =' endif goto 6040 endif 5001 format(14x,a) C C---rm() d-parameter C if(nparf(n,1).eq.-5)then I=nparf(n,2) if(i.le.3)then WRITE(NAMPAR(N),5001)'d_'//parnm(i)//' =' else if(i.eq.7) WRITE(NAMPAR(N),'(8x,a)')'d_a = d_b =' if(i.eq.8) WRITE(NAMPAR(N),'(8x,a)')'d_a = d_c =' if(i.eq.9) WRITE(NAMPAR(N),'(8x,a)')'d_b = d_c =' if(i.eq.10)WRITE(NAMPAR(N),'(2x,a)')'d_a = d_b = d_c =' endif goto 6040 endif C C...Laurie hydrogen atom type correction C if(nparf(n,1).eq.-6)then I=nparf(n,2) j=ncart(i,1) jj=ncart(i,2) WRITE(NAMPAR(N),6055)j,jj 6055 FORMAT( ' Laurie d_H(',I2,',',I2,') =') iok=0 if(cartrf(j,4).gt.1.d0.and.cartrf(j,4).lt.3.02d0)iok=1 if(cartrf(jj,4).gt.1.d0.and.cartrf(jj,4).lt.3.02d0)iok=1 if(iok.eq.0)then WRITE(NAMPAR(N),508)j,jj 508 FORMAT( ' Laurie d_R(',I2,',',I2,') =') write(*,500)j,jj write(3,500)j,jj endif goto 6040 endif 500 format(' ***** WARNING------> Laurie parameter for bond (', * i2,',',i2,'): no hydrogens involved') C C...Laurie hydrogen atom type correction defined in terms of projection C angle (only for negative values of delta_H) C if(nparf(n,1).eq.-7)then I=nparf(n,2) j=ncart(i,1) jj=ncart(i,2) WRITE(NAMPAR(N),6075)j,jj 6075 FORMAT( ' Laurie A.r(',I2,',',I2,') =') iok=0 if(cartrf(j,4).gt.1.d0.and.cartrf(j,4).lt.3.02d0)iok=1 if(cartrf(jj,4).gt.1.d0.and.cartrf(jj,4).lt.3.02d0)iok=1 if(iok.eq.0)then WRITE(NAMPAR(N),578)j,jj 578 FORMAT( ' Laurie A.R(',I2,',',I2,') =') write(*,500)j,jj write(3,500)j,jj endif goto 6040 endif C C...Laurie type correction to a bond angle C if(nparf(n,1).eq.-9)then I=nparf(n,2) j=ncart(i,1) jj=ncart(i,2) jjj=ncart(i,3) WRITE(NAMPAR(N),6065)j,jj,jjj 6065 FORMAT( ' Lt d_HA(',I2,',',I2,',',I2,') =') goto 6040 endif C C...deuteration shrinkage C if(nparf(n,1).eq.-8)then I=nparf(n,2) j=ncart(i,2) WRITE(NAMPAR(N),6056)i,j 6056 FORMAT( ' delta_D(',I2,',',I2,') =') goto 6040 endif C C...structural internal C if(value.eq.0.d0)VPARF(N)=CARTDF(NPARF(N,1),NPARF(N,2)) I=NPARF(N,1) if(nparf(n,2).eq.1) * WRITE(NAMPAR(N),6035)NCART(I,1),NCART(I,2) if(nparf(n,2).eq.2) * WRITE(NAMPAR(N),6036)NCART(I,1),NCART(I,2),NCART(I,3) if(nparf(n,2).eq.3) * WRITE(NAMPAR(N),6037)NCART(I,1),NCART(I,2),NCART(I,3), * NCART(I,4) 6035 FORMAT(10x, 'R(',I2,',',I2,') =') 6036 FORMAT( 7x, 'A(',I2,',',I2,',',I2,') =') 6037 FORMAT( 4x,'D(',I2,',',I2,',',I2,',',I2,') =') C 6040 if(ifit(n).eq.1)WRITE(3,6041)NAMPAR(N),VPARF(N) 6041 FORMAT(1X,A,F13.6) C C...read any repetitions C IF(NPARF(N,3).NE.0)THEN READ(2,'(a)',err=6009)LINE nlin=nlin+1 READ(LINE,6047,err=6008)TEXT,(NREP(N,NN),NN=1,NPARF(N,3)) if(ifit(n).eq.1) * WRITE(3,6046)NPARF(N,3),(NREP(N,NN),NN=1,NPARF(N,3)) 6046 FORMAT(27X,', and at',I2,' more atom(s): ',20I3) c c...check for errors if(nparf(n,1).eq.-6.or.nparf(n,1).eq.-7)then do 503 nn=1,nparf(n,3) I=nrep(n,nn) j=ncart(i,1) jj=ncart(i,2) c if(cartrf(j,4).lt.1.d0.or.cartrf(jj,4).lt.1.d0)goto 505 iok=0 if(cartrf(j,4).gt.1.d0.and.cartrf(j,4).lt.3.02d0)iok=1 if(cartrf(jj,4).gt.1.d0.and.cartrf(jj,4).lt.3.02d0)iok=1 if(iok.eq.0)then write(*,500)j,jj write(3,500)j,jj endif 503 continue endif c ENDIF 6047 FORMAT(A,20I2) 6004 CONTINUE C IF(NFIX.GT.0)THEN WRITE(3,'(1x/1x,a/14x,17(1h-))') * ' Fixed parameters: ' do 6070 nn=1,npar if(ifit(nn).eq.0)then WRITE(3,6041)NAMPAR(NN),VPARF(NN) IF(NPARF(NN,3).NE.0) * WRITE(3,6046)NPARF(NN,3),(NREP(NN,NNN),NNN=1,NPARF(NN,3)) endif 6070 continue ENDIF C C - - - - - - - - - - - - - C C...Read spectroscopic constants: nisot,nconst,Bobs,deltaI,deltaB C 6701 READ(2,'(a)',err=6009)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE goto 6701 endif C NISOT=1 READ(LINE,6003,ERR=6008)TEXT,NOBSC WRITE(3,6011)' NO OF CONSTANTS TO FIT TO: ',NOBSC DO 6005 N=1,NOBSC READ(2,'(A)',ERR=6009)LINE nlin=nlin+1 READ(LINE,6003,ERR=6008)TEXT,(NCOBS(N,I),I=1,2),VCOBS(N,3), * DELTAI,DELTAB c if(DELTAB.ne.0.d0.and.DELTAI.ne.0.d0)then write(*,'(1x//'' **** ERROR: only one of dI and dB '', * ''can be non-zero''//18x,2i2,3f12.5//)') * (NCOBS(N,I),I=1,2),VCOBS(N,3),DELTAI,DELTAB stop endif if(DELTAB.ne.0.d0.and.NCOBS(N,1).gt.3)then write(*,'(1x//'' **** ERROR: use of dB only allowed '', * ''for A,B,C. Correct this line:''//a//)')line stop endif c VCOBS(N,1)=CONV/VCOBS(N,3) if(ncobs(n,1).eq.4)vcobs(n,1)=4.d0*vcobs(n,1) if(ncobs(n,1).eq.5)vcobs(n,1)=4.d0*vcobs(n,1) if(ncobs(n,1).ge.6)then vcobs(n,1)=vcobs(n,3) vcobs(n,3)=conv/vcobs(n,1) endif C IF(NCOBS(N,2).LE.NISOT)GOTO 6030 C C...entry for a new isotopic species C NISOT=NCOBS(N,2) c IF(DELTAI.NE.0.D0)THEN VCOBS(N,1)=VCOBS(N,1)+DELTAI WRITE(3,7038)NISOT,CNAME(NCOBS(N,1)), * VCOBS(N,3),CONV/VCOBS(N,1),DELTAI VCOBS(N,3)=CONV/VCOBS(N,1) 7038 FORMAT(' Isotopic species',I3,': ',A,'exp =',F12.5,3X, * 'corr =',F12.5,' (dI =',F12.5,')') ENDIF c IF(DELTAB.NE.0.0d0)THEN WRITE(3,7138)NISOT,CNAME(NCOBS(N,1)), * VCOBS(N,3),VCOBS(N,3)+DELTAB,DELTAB VCOBS(N,3)=VCOBS(N,3)+DELTAB VCOBS(N,1)=CONV/VCOBS(N,3) 7138 FORMAT(' Isotopic species',I3,': ',A,'exp =',F12.5,3X, * 'corr =',F12.5,' (dB =',F12.5,')') ENDIF c IF(DELTAI.eq.0.d0.AND.DELTAB.eq.0.0d0)THEN WRITE(3,6038)NISOT,CNAME(NCOBS(N,1)),VCOBS(N,3) 6038 FORMAT(' Isotopic species',I3,': ',A,' =',F12.5) ENDIF GOTO 6005 c c...entry for an existing isotopic species c 6030 IF(DELTAI.NE.0.0D0)THEN VCOBS(N,1)=VCOBS(N,1)+DELTAI WRITE(3,7031)CNAME(NCOBS(N,1)), * VCOBS(N,3),CONV/VCOBS(N,1),DELTAI if(vcobs(n,1).ne.0.d0)then VCOBS(N,3)=CONV/VCOBS(N,1) else vcobs(n,3)=1.d+15 endif 7031 FORMAT(24X,A,'exp =',F12.5,3X, * 'corr =',F12.5,' (dI =',F12.5,')') ENDIF c IF(DELTAB.NE.0.0d0)THEN WRITE(3,7131)CNAME(NCOBS(N,1)), * VCOBS(N,3),VCOBS(N,3)+DELTAB,DELTAB VCOBS(N,3)=VCOBS(N,3)+DELTAB VCOBS(N,1)=CONV/VCOBS(N,3) 7131 FORMAT(24X,A,'exp =',F12.5,3X, * 'corr =',F12.5,' (dB =',F12.5,')') ENDIF c IF(DELTAI.eq.0.d0.AND.DELTAB.eq.0.0d0)THEN WRITE(3,6031)CNAME(NCOBS(N,1)),VCOBS(N,3) 6031 FORMAT(24X,A,' =',F12.5) ELSE ENDIF c 6005 CONTINUE C C - - - - - - - - - - - - - C C...Read definitions of substituted species C NUMCHG(1)=0 IF(NISOT.EQ.1)GOTO 6006 write(3,'(1x)') DO 6007 N=2,NISOT 6704 READ(2,'(a)',err=6009)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE goto 6704 endif READ(LINE,6003,ERR=6008)TEXT,NUMCHG(N) WRITE(3,6042)N 6042 FORMAT(' ISOTOPIC SPECIES',I3,', changes from ', 1 'parent species:') DO 6017 NN=1,NUMCHG(N) READ(2,'(A)',ERR=6009)LINE nlin=nlin+1 READ(LINE,6093,ERR=6008) * TEXT,(NCHANG(N,NN,I),I=1,2),VCHANG(N,NN) WRITE(3,6043)TEXT,(NCHANG(N,NN,I),I=1,2),VCHANG(N,NN) 6043 FORMAT(1X,A,2I3,F17.7) 6017 CONTINUE 6007 CONTINUE 6093 FORMAT(A,2I2,F20.7) GOTO 6006 C 6008 WRITE(*,'(1X//'' **** ERROR on reading from line'',i5, * '' in the input file:''//1x,A,A//)')nlin,LINE,CHAR(7) STOP c 6009 WRITE(*,'(1X//'' **** ERROR on reading line'',i5, * '' of the input file:''//1x,A,A//)')nlin+1,FILNAM,CHAR(7) STOP C 6006 CLOSE(2) c RETURN END C C_____________________________________________________________________________ C SUBROUTINE LMUSER(X,A,Y,DYDA,NA) C C User supplied routine for the Levenberg-Marquardt nonlinear least-squares C fit. For x=X and using the NA constants in A this returns the C calculated value of y in Y and partial derivatives of y wrt all the C constants in DYDA. Since the necessary calculations have already been C performed by SETUP only a transfer of data between various arrays is C required here. C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=50,maxiso=30,maxcon=3*maxiso,maxrep=20, * maxchg=10,maxato=300) C COMMON /DERIV/DER(MAXCON,maxpar) COMMON /INPUT/CARTRF(MAXATO,4),NPARF(maxpar,3), * NREP(maxpar,maxrep),VPARF(maxpar),NCOBS(MAXCON,2), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2), * VCHANG(maxiso,maxchg),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,COMENT REAL*8 X,A(NA),Y,DYDA(NA) CHARACTER*20 NAMPAR(maxpar) CHARACTER*78 COMENT C N=NINT(X) DO 1 I=1,NA DYDA(I)=DER(N,I) 1 CONTINUE Y=VCOBS(N,2) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SETUP(A,MA) C C This routine provides an interface to the Levenberg-Marquardt set of C routines - it evaluates calculated values and derivatives (wrt all C parameters) for all constants. In this way it only remains for the user C supplied routine FUNCS in the ML scheme to extract the appropriate C equation for a given constant. C C A - vector of length MA containing the current constants to be used for C evaluations of Y.calc (MFIT=MA of which are to be fitted) C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=50,maxiso=30,maxcon=3*maxiso,maxrep=20, * maxchg=10,maxato=300) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /INPUT/CARTRF(MAXATO,4),NPARF(maxpar,3), * NREP(maxpar,maxrep),VPARF(maxpar),NCOBS(MAXCON,2), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2), * VCHANG(maxiso,maxchg),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,COMENT COMMON /DERIV/DER(MAXCON,maxpar) REAL*8 A(MA) CHARACTER*20 NAMPAR(maxpar) CHARACTER*78 COMENT C C...zeroth run C izero=0 CALL SETUP1(izero,izero,A,MA) <------- C C...derivative runs for all parameters of fit in turn C N is position in A etc. C NV is position in VPARF etc. C n=0 DO 1 NV=1,NPAR if(ifit(nv).eq.0)goto 1 n=n+1 CALL SETUP1(N,nv,A,MA) <------- 1 CONTINUE C WRITE(3,4) 4 FORMAT(1X/'I.obs, I.obs-calc and derivatives:'/) DO 2 N=1,NOBSC WRITE(3,3)VCOBS(N,1),VCOBS(N,1)-VCOBS(N,2), 1 (DER(N,NN),NN=1,MA) 2 CONTINUE 3 FORMAT(1X,2F10.5,' ...',6F9.3/(25X,6F9.3)) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SETUP1(NCALL,NV,A,MA) C C NCALL=0 calculate values for all measured spectroscopic constants C NCALL=i calculate derivatives of all spectroscopic constants with respect C to the i'th parameter of fit C NCALL=-1 calculate values for all measured spectroscopic constants but C only from the direct structural internal and without using C any of those parameters which contribute to the eps.alpha C contribution to I.alpha C C A - vector of length MA containing only the parameters to be fitted, C (MFIT=MA) C NV - position of the current parameter of fit as defined by NCALL among C all of the declared parameters C VPARF - contains values of all the constants to be used for evaluations C of Y.calc - these are updated in M/P from the latest values in A C C NOTE - indexing of A,AOLD,ERRPAR,COVAR is 1..MA (MA.le.NPAR) C indexing of VPARF,NPARF,NAMPAR is 1..NPAR C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=50,maxiso=30,maxcon=3*maxiso,maxrep=20, * maxchg=10,maxato=300) parameter (cang=57.2957795131d0) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /ANCIL/RMC(6),RMD(3),rmci(6),rmdi(3) COMMON /coords/CORDP(MAXATO,3) COMMON /INPUT/CARTRF(MAXATO,4),NPARF(maxpar,3), * NREP(maxpar,maxrep),VPARF(maxpar),NCOBS(MAXCON,2), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2), * VCHANG(maxiso,maxchg),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,COMENT COMMON /DERIV/DER(MAXCON,maxpar) COMMON /EPSIL/ristr(maxcon,3),ritot(maxcon,3) REAL*8 DELTA(9) REAL*8 A(MA),rlaur(maxpar),rlaurr(maxpar,maxrep) CHARACTER*20 NAMPAR(maxpar) CHARACTER*78 COMENT DATA DELTA/0.00001D0,0.001D0,0.001D0, r,a,t * 0.00001D0,0.00001D0, c,d * 0.00001D0,0.01000D0, dh,dhpra * 0.00001D0,0.01000D0/ sr,dlang C C C...Main loop over all isotopic species, NC counts the number of C spectroscopic constants that have been dealt with C NC=1 DO 1 N=1,NISOT C C...set up the current CART data for calculation by transferring into CARTDF: C 1/ the reference configuration held in CARTRF C 2/ the current parameters of fit in A C C...Transfer the reference configuration C DO 3 NN=1,NATMS DO 3 NNN=1,4 CARTDF(NN,NNN)=CARTRF(NN,NNN) 3 CONTINUE C C...Replace reference internals by current parameters of fit C kk=0 DO 20 NN=1,NPAR if(nparf(nn,1).lt.0)goto 20 if(ifit(nn).eq.0)then CARTDF(NPARF(NN,1),NPARF(NN,2))=VPARF(NN) IF(NPARF(NN,3).NE.0)THEN DO 22 NNN=1,NPARF(NN,3) CARTDF(iabs(NREP(NN,NNN)),NPARF(NN,2))=VPARF(NN) if(nrep(nn,nnn).lt.0) * CARTDF(iabs(NREP(NN,NNN)),NPARF(NN,2))=-VPARF(NN) 22 CONTINUE ENDIF else kk=kk+1 CARTDF(NPARF(NN,1),NPARF(NN,2))=A(kk) IF(NPARF(NN,3).NE.0)THEN DO 32 NNN=1,NPARF(NN,3) CARTDF(iabs(NREP(NN,NNN)),NPARF(NN,2))=A(kk) if(nrep(nn,nnn).lt.0) * CARTDF(iabs(NREP(NN,NNN)),NPARF(NN,2))=-A(kk) 32 CONTINUE ENDIF endif 20 CONTINUE C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...Set up the current RM(1) and RM(2) coefficients in RMC and RMD C from any such constants in A C do 101 nn=1,3 rmc(nn)= 0.d0 rmc(nn+3)=0.d0 rmd(nn)= 0.d0 101 continue if(ncall.lt.0)goto 415 c kk=0 do 100 nn=1,npar if(ifit(nn).eq.1)kk=kk+1 if(nparf(nn,1).eq.-4)then i=nparf(nn,2) c if(i.le.6)then if(ifit(nn).eq.1)then rmc(i)=a(kk) else rmc(i)=vparf(nn) endif endif c if(i.eq.7)then if(ifit(nn).eq.1)then rmc(1)=a(kk) rmc(2)=a(kk) else rmc(1)=vparf(nn) rmc(2)=vparf(nn) endif endif c if(i.eq.8)then if(ifit(nn).eq.1)then rmc(1)=a(kk) rmc(3)=a(kk) else rmc(1)=vparf(nn) rmc(3)=vparf(nn) endif endif c if(i.eq.9)then if(ifit(nn).eq.1)then rmc(2)=a(kk) rmc(3)=a(kk) else rmc(2)=vparf(nn) rmc(3)=vparf(nn) endif endif c if(i.eq.10)then if(ifit(nn).eq.1)then rmc(1)=a(kk) rmc(2)=a(kk) rmc(3)=a(kk) else rmc(1)=vparf(nn) rmc(2)=vparf(nn) rmc(3)=vparf(nn) endif endif c endif C if(nparf(nn,1).eq.-5)then i=nparf(nn,2) c if(i.le.6)then if(ifit(nn).eq.1)then rmd(i)=a(kk) else rmd(i)=vparf(nn) endif endif c if(i.eq.7)then if(ifit(nn).eq.1)then rmd(1)=a(kk) rmd(2)=a(kk) else rmd(1)=vparf(nn) rmd(2)=vparf(nn) endif endif c if(i.eq.8)then if(ifit(nn).eq.1)then rmd(1)=a(kk) rmd(3)=a(kk) else rmd(1)=vparf(nn) rmd(3)=vparf(nn) endif endif c if(i.eq.9)then if(ifit(nn).eq.1)then rmd(2)=a(kk) rmd(3)=a(kk) else rmd(2)=vparf(nn) rmd(3)=vparf(nn) endif endif c if(i.eq.10)then if(ifit(nn).eq.1)then rmd(1)=a(kk) rmd(2)=a(kk) rmd(3)=a(kk) else rmd(1)=vparf(nn) rmd(2)=vparf(nn) rmd(3)=vparf(nn) endif endif c endif 100 continue C 415 IF(N.EQ.1)GOTO 414 C C...implement changes for the isotopic species C DO 4 NN=1,NUMCHG(N) I=NCHANG(N,NN,1) J=NCHANG(N,NN,2) IF(J.NE.4)GOTO 8 CARTDF(I,J)=VCHANG(N,NN) GOTO 4 8 CARTDF(I,J)=CARTDF(I,J)+VCHANG(N,NN) 4 CONTINUE C C...evaluate current molecular mass C 414 if(ncall.eq.-1)goto 2 ww= 0.d0 wwp=0.d0 do 201 nnn=1,natms ww= ww +cartdf(nnn,4) wwp=wwp+cartrf(nnn,4) 201 continue C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...add Laurie correction if applicable C kk=0 DO 200 nn=1,npar if(ifit(nn).eq.1)kk=kk+1 c c...standard linear correction c if(nparf(nn,1).eq.-6)then c c...correction to first bond i=nparf(nn,2) rmh=cartdf(i,4) rmha=rmh j=ncart(i,2) if(rmh.gt.3.017d0.or.rmh.lt.1.0d0)then rmh=cartdf(j,4) if(rmh.gt.3.017d0)rmh=0.5d0*(rmh+rmha) endif if(j.gt.i)i=j if(ifit(nn).eq.1)then rlaur(nn)=a(kk)*dsqrt(ww/(rmh*(ww-rmh))) else rlaur(nn)=vparf(nn)*dsqrt(ww/(rmh*(ww-rmh))) endif cartdf(i,1)=cartdf(i,1)+rlaur(nn) C C...deal with repetitions if(nparf(nn,3).ne.0)then do 215 nnn=1,nparf(nn,3) i=nrep(nn,nnn) rmh=cartdf(i,4) rmha=rmh j=ncart(i,2) if(rmh.gt.3.017d0.or.rmh.lt.1.d0)then rmh=cartdf(j,4) if(rmh.gt.3.017d0)rmh=0.5d0*(rmh+rmha) endif if(j.gt.i)i=j if(ifit(nn).eq.1)then rlaurr(nn,nnn)=a(kk)*dsqrt(ww/(rmh*(ww-rmh))) else rlaurr(nn,nnn)=vparf(nn)*dsqrt(ww/(rmh*(ww-rmh))) endif cartdf(i,1)=cartdf(i,1)+rlaurr(nn,nnn) 215 continue endif endif C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c...Laurie correction defined in terms of an angle c if(nparf(nn,1).eq.-7)then c c...apply correction to first bond i=nparf(nn,2) rmh =cartdf(i,4) rmhp=cartrf(i,4) rmha =rmh rmhap=rmhp j=ncart(i,2) if(rmh.gt.3.017d0.or.rmh.lt.1.0d0)then rmh =cartdf(j,4) rmhp=cartrf(j,4) if(rmh.gt.3.017d0)then rmh =0.50d0*(rmh +rmha) rmhp=0.50d0*(rmhp+rmhap) c rmh =1.00d0*(rmh +rmha) c rmhp=1.00d0*(rmhp+rmhap) endif endif if(j.gt.i)i=j if(ifit(nn).eq.1)then costh=dcos(a(kk)/cang) rlaur(nn)=(costh-1.d0)*cartdf(i,1) * *dsqrt(rmhp*(wwp-rmhp)/wwp) * *dsqrt( ww/(rmh*(ww-rmh) )) else costh=dcos(vparf(nn)/cang) rlaur(nn)=(costh-1.d0)*cartdf(i,1) * *dsqrt(rmhp*(wwp-rmhp)/wwp) * *dsqrt( ww/(rmh*(ww-rmh) )) endif cosths=(cartdf(i,1)+rlaur(nn))/cartdf(i,1) cartdf(i,1)=cartdf(i,1)+rlaur(nn) C C...deal with repetitions if(nparf(nn,3).ne.0)then do 419 nnn=1,nparf(nn,3) i=nrep(nn,nnn) j=ncart(i,2) if(j.gt.i)i=j rlaurr(nn,nnn)=(cosths-1.d0)*cartdf(i,1) cartdf(i,1)=cartdf(i,1)+rlaurr(nn,nnn) 419 continue endif endif 200 continue C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...Laurie type correction for an angle (do not use - only for experimenting with C new terms and not completed!) C kk=0 DO 220 nn=1,npar if(ifit(nn).eq.1)kk=kk+1 if(nparf(nn,1).eq.-9)then c c...correction to first angle i=nparf(nn,2) rmh=cartdf(i,4) if(ifit(nn).eq.1)then rlaur(nn)=a(kk)*dsqrt(ww/(rmh*(ww-rmh))) else rlaur(nn)=vparf(nn)*dsqrt(ww/(rmh*(ww-rmh))) endif cartdf(i,2)=cartdf(i,2)+rlaur(nn) C C...deal with repetitions if(nparf(nn,3).ne.0)then do 221 nnn=1,nparf(nn,3) i=nrep(nn,nnn) rmh=cartdf(i,4) if(ifit(nn).eq.1)then rlaur(nn)=a(kk)*dsqrt(ww/(rmh*(ww-rmh))) else rlaur(nn)=vparf(nn)*dsqrt(ww/(rmh*(ww-rmh))) endif cartdf(i,2)=cartdf(i,2)+rlaur(nn) 221 continue endif endif 220 continue C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...add the deuteration stretch if the relevant bond is substituted C kk=0 do 235 nn=1,npar if(ifit(nn).eq.1)kk=kk+1 if(nparf(nn,1).eq.-8)then i=nparf(nn,2) j=ncart(i,2) if( (cartdf(i,4).gt.1.d0).and. * (cartdf(i,4).lt.3.03d0) .and. * (cartdf(i,4)-cartrf(i,4).ne.0.d0) ) * then if(ifit(nn).eq.1)then cartdf(i,1)=cartdf(i,1)+a(kk) else cartdf(i,1)=cartdf(i,1)+vparf(nn) endif if(idump.eq.1.and.ncall.eq.0)write(3,251)n,i,j 251 format(30x,'Isotope',I3,', bond stretch for bond:',2i3) endif if( (cartdf(j,4).gt.1.d0).and. * (cartdf(j,4).lt.3.03d0).and. * (cartdf(j,4)-cartrf(j,4).ne.0.d0) ) * then if(ifit(nn).eq.1)then cartdf(i,1)=cartdf(i,1)+a(kk) else cartdf(i,1)=cartdf(i,1)+vparf(nn) endif if(idump.eq.1.and.ncall.eq.0)write(3,251)n,i,j endif C C...deal with repetitions IF(NPARF(NN,3).NE.0)THEN DO 236 NNN=1,NPARF(NN,3) i=nrep(nn,nnn) j=ncart(i,2) if( (cartdf(i,4).gt.1.d0).and. * (cartdf(i,4).lt.3.03d0) .and. * (cartdf(i,4)-cartrf(i,4).ne.0.d0) ) * then if(ifit(nn).eq.1)then cartdf(i,1)=cartdf(i,1)+a(kk) else cartdf(i,1)=cartdf(i,1)+vparf(nn) endif if(idump.eq.1.and.ncall.eq.0)write(3,251)n,i,j endif if( (cartdf(j,4).gt.1.d0).and. * (cartdf(j,4).lt.3.03d0) .and. * (cartdf(j,4)-cartrf(j,4).ne.0.d0) ) * then if(ifit(nn).eq.1)then cartdf(i,1)=cartdf(i,1)+a(kk) else cartdf(i,1)=cartdf(i,1)+vparf(nn) endif if(idump.eq.1.and.ncall.eq.0)write(3,251)n,i,j endif 236 CONTINUE ENDIF endif 235 continue C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C 14 IF(NCALL.LE.0)GOTO 2 C C...increment the current fitted parameter for derivative calculation C this is parameter NV in list of all parameters C NCAll in list of parameters of fit C I=NPARF(NV,1) J=NPARF(NV,2) C C...gradient for c's C if(i.eq.-4)then if(j.le.6)rmc(j)=rmc(j)+delta(4) if(j.eq.7)then rmc(1)=rmc(1)+delta(4) rmc(2)=rmc(2)+delta(4) endif if(j.eq.8)then rmc(1)=rmc(1)+delta(4) rmc(3)=rmc(3)+delta(4) endif if(j.eq.9)then rmc(2)=rmc(2)+delta(4) rmc(3)=rmc(3)+delta(4) endif if(j.eq.10)then rmc(1)=rmc(1)+delta(4) rmc(2)=rmc(2)+delta(4) rmc(3)=rmc(3)+delta(4) endif goto 2 endif C C...gradient for d's C if(i.eq.-5)then if(j.le.3)rmd(j)=rmd(j)+delta(5) if(j.eq.7)then rmd(1)=rmd(1)+delta(5) rmd(2)=rmd(2)+delta(5) endif if(j.eq.8)then rmd(1)=rmd(1)+delta(5) rmd(3)=rmd(3)+delta(5) endif if(j.eq.9)then rmd(2)=rmd(2)+delta(5) rmd(3)=rmd(3)+delta(5) endif if(j.eq.10)then rmd(1)=rmd(1)+delta(5) rmd(2)=rmd(2)+delta(5) rmd(3)=rmd(3)+delta(5) endif goto 2 endif C C...gradient for standard Laurie contribution C if(i.eq.-6)then i=nparf(NV,2) rmh=cartdf(i,4) rmha=rmh j=ncart(i,2) if(rmh.gt.3.017d0.or.rmh.lt.1.d0)then rmh=cartdf(j,4) if(rmh.gt.3.017d0)rmh=0.5d0*(rmh+rmha) endif if(j.gt.i)i=j cartdf(i,1)=cartdf(i,1)+delta(6)*dsqrt(ww/(rmh*(ww-rmh))) C if(nparf(nv,3).ne.0)then do 216 nnn=1,nparf(nv,3) i=nrep(nv,nnn) rmh=cartdf(i,4) rmha=rmh j=ncart(i,2) if(rmh.gt.3.017d0.or.rmh.lt.1.d0)then rmh=cartdf(j,4) if(rmh.gt.3.017d0)rmh=0.5d0*(rmh+rmha) endif if(j.gt.i)i=j cartdf(i,1)=cartdf(i,1)+delta(6)*dsqrt(ww/(rmh*(ww-rmh))) 216 continue endif goto 2 endif C C...gradient for Laurie contribution defined by projection angle C if(i.eq.-7)then i=nparf(NV,2) rmh =cartdf(i,4) rmhp=cartrf(i,4) rmha =rmh rmhap=rmhp j=ncart(i,2) if(rmh.gt.3.017d0.or.rmh.lt.1.d0)then rmh =cartdf(j,4) rmhp=cartrf(j,4) if(rmh.gt.3.017d0)then rmh= 0.50d0*(rmh +rmha) rmhp=0.50d0*(rmhp+rmhap) endif endif if(j.gt.i)i=j cartdf(i,1)=cartdf(i,1)-rlaur(nv) costhd=dcos( (a(ncall)+delta(7)) /cang) rlarrr=(costhd-1.d0)*cartdf(i,1) * *dsqrt(rmhp*(wwp-rmhp)/wwp) * *dsqrt( ww/(rmh*(ww-rmh) )) cosths=(cartdf(i,1)+rlarrr)/cartdf(i,1) cartdf(i,1)=cartdf(i,1)+rlarrr C if(nparf(nv,3).ne.0)then do 416 nnn=1,nparf(nv,3) i=nrep(nv,nnn) j=ncart(i,2) if(j.gt.i)i=j cartdf(i,1)=cartdf(i,1)-rlaurr(nv,nnn) rlarrr=(cosths-1.d0)*cartdf(i,1) cartdf(i,1)=cartdf(i,1)+rlarrr 416 continue endif goto 2 endif C C...gradient for Laurie type contribution to bond angle C if(i.eq.-9)then i=nparf(NV,2) j=ncart(i,3) rmh=cartdf(i,4) cartdf(i,2)=cartdf(i,2)+delta(7)*dsqrt(ww/(rmh*(ww-rmh))) C if(nparf(nv,3).ne.0)then do 226 nnn=1,nparf(nv,3) i=nrep(nv,nnn) rmh=cartdf(i,4) cartdf(i,2)=cartdf(i,2)+delta(7)*dsqrt(ww/(rmh*(ww-rmh))) 226 continue endif goto 2 endif C C...gradient for deuteration stretch C if(i.eq.-8)then i=nparf(nv,2) j=ncart(i,2) if( (cartdf(i,4).gt.1.d0).and. * (cartdf(i,4).lt.3.03d0) .and. * (cartdf(i,4)-cartrf(i,4).ne.0.d0) ) * cartdf(i,1)=cartdf(i,1)+delta(8) if( (cartdf(j,4).gt.1.d0).and. * (cartdf(j,4).lt.3.03d0) .and. * (cartdf(j,4)-cartrf(j,4).ne.0.d0) ) * cartdf(i,1)=cartdf(i,1)+delta(8) C IF(NPARF(nv,3).NE.0)THEN DO 240 NNN=1,NPARF(NV,3) i=nrep(nv,nnn) j=ncart(i,2) if( (cartdf(i,4).gt.1.d0).and. * (cartdf(i,4).lt.3.03d0) .and. * (cartdf(i,4)-cartrf(i,4).ne.0.d0) ) * cartdf(i,1)=cartdf(i,1)+delta(8) if( (cartdf(j,4).gt.1.d0).and. * (cartdf(j,4).lt.3.03d0) .and. * (cartdf(j,4)-cartrf(j,4).ne.0.d0) ) * cartdf(i,1)=cartdf(i,1)+delta(8) 240 CONTINUE ENDIF goto 2 endif C C...gradient for a standard internal: bond, bond angle, dihedral angle C CARTDF(I,J)=CARTDF(I,J)+DELTA(J) IF(NPARF(NV,3).NE.0)THEN DO 23 NNN=1,NPARF(NV,3) I=iabs(NREP(NV,NNN)) CARTDF(I,J)=CARTDF(I,J)+DELTA(J) if(nrep(nv,nnn).lt.0)cartdf(i,j)=cartdf(i,j)-2*delta(j) 23 CONTINUE ENDIF C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...calculate moments of inertia for this (=Nth) isotopomer C 2 IF(N.EQ.1.and.ncall.eq.0)then niter=niter+1 WRITE(3,'(1X/1x,78(1h-)/'' ITERATION'',i5/1x,78(1h-)/)') * niter endif C CALL MOMIN(n,ncall) <------- C if(ncall.eq.-1)then ristr(n,1)=pmi(1,1) ristr(n,2)=pmi(2,1) ristr(n,3)=pmi(3,1) endif if(ncall.eq.0)then ritot(n,1)=pmi(1,1) ritot(n,2)=pmi(2,1) ritot(n,3)=pmi(3,1) endif C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...write key results if this is the zeroth run C IF(NCALL.eq.0)THEN if(pmi(1,2).lt.1.d6)then WRITE(3,9)N,(PMI(NN,2),NN=1,3),(PMI(nn,1),nn=1,3) else WRITE(3,99)N,(PMI(NN,2),NN=2,3),(PMI(nn,1),nn=1,3) endif c IF(.NOT. * (RMCI(1).EQ.0.D0.AND.RMCI(2).EQ.0.D0.AND.RMCI(3).EQ.0.D0)) * WRITE(3,91)(rmci(nn),nn=1,3) IF(.NOT. * (RMCI(4).EQ.0.D0.AND.RMCI(5).EQ.0.D0.AND.RMCI(6).EQ.0.D0)) * write(3,92)(rmci(nn),nn=4,6) IF(.NOT. * (RMDI(1).EQ.0.D0.AND.RMDI(2).EQ.0.D0.AND.RMDI(3).EQ.0.D0)) * WRITE(3,93)(rmdi(nn),nn=1,3) c 9 FORMAT(' Isot.species',I3,': A =',F12.5,5X,'B =', * F12.5,5X,'C =',F12.5/ * 21x,'I.a =',F12.5,3x,'I.b =',F12.5,3x,'I.c =',F12.5) 99 FORMAT(' Isot.species',I3,': A =',12(1H*),5X,'B =', * F12.5,5X,'C =',F12.5/ * 21x,'I.a =',F12.5,3x,'I.b =',F12.5,3x,'I.c =',F12.5) 91 format( * 10x,'incr. from c.a =',F12.5,3x,'c.b =',F12.5,3x,'c.c =',F12.5) 92 FORMAT( * 10x,'incr. from c.ab=',F12.5,3x,'c.ac=',F12.5,3x,'c.bc=',F12.5) 93 format( * 10x,'incr. from d.a =',F12.5,3x,'d.b =',F12.5,3x,'d.c =',F12.5) c do 280 nn=1,npar if(nparf(nn,1).eq.-6.or.nparf(nn,1).eq.-7)then i=nparf(nn,2) if(rlaur(nn).eq.0.d0)goto 280 if(nparf(nn,1).eq.-6)then write(3,94)i,ncart(i,2),cartdf(i,1)-rlaur(nn), * rlaur(nn),cartdf(i,1) else angle=cang*dacos(cartdf(i,1)/(cartdf(i,1)-rlaur(nn))) write(3,95)i,ncart(i,2),cartdf(i,1)-rlaur(nn), * angle,cartdf(i,1) endif if(nparf(nn,3).ne.0)then do 217 nnn=1,nparf(nn,3) i=nrep(nn,nnn) if(nparf(nn,1).eq.-6)then write(3,94)i,ncart(i,2),cartdf(i,1)-rlaurr(nn,nnn), * rlaurr(nn,nnn),cartdf(i,1) else angle=cang* * dacos( cartdf(i,1)/(cartdf(i,1)-rlaurr(nn,nnn)) ) write(3,95)i,ncart(i,2), * cartdf(i,1)-rlaurr(nn,nnn), * angle,cartdf(i,1) endif 217 continue endif endif 280 continue 94 FORMAT(12x,' r(',i2,',',i2,') =',F12.5,'(rigid) +', * F11.5,'(Laurie) =',F10.5) 95 FORMAT(12x,' r(',i2,',',i2,') =',F12.5,'(rigid) * cos(', * F10.5,') =',F10.5) ENDIF C C...extract the required constant or derivative C 5 I=NCOBS(NC,1) J=NCOBS(NC,2) IF(J.GT.N)GOTO 1 C IF(NCALL.eq.0)then VCOBS(NC,2)=PMI(I,1) VCOBS(NC,4)=PMI(I,2) ELSE ii=nparf(nv,1) if(ii.gt.0)then DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/DELTA(NPARF(NV,2)) else if(ii.eq.-4)DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/delta(4) if(ii.eq.-5)DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/delta(5) if(ii.eq.-6)DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/delta(6) if(ii.eq.-7)DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/delta(7) if(ii.eq.-8)DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/delta(8) if(ii.eq.-9)DER(NC,NCALL)=(PMI(I,1)-VCOBS(NC,2))/delta(9) endif ENDIF C 16 NC=NC+1 IF(NC.GT.NOBSC)GOTO 17 GOTO 5 C 1 CONTINUE C C 17 RETURN END C C_____________________________________________________________________________ C SUBROUTINE MOMIN(isotop,ncall) C C Principal moments of inertia of a molecule from its CART type C set of definitions C C ISOTOP declares the number of the isotopic species for this C evaluation C (which may be for the zeroth run or for a derivative run for C any of the fitted constants) C NCALL specifies that MOMIN has been called during: C = 0 the zeroth run C >1 a derivative run C -1 first step of eps evaluation C C MOMIN returns, for i=1,2,3,4,5,6,7,8 for a,b,c,b+c,a+b,Pa,Pb,Pc: C C PMI(i,1) - moments of inertia C PMI(i,2) - rotational constants C RMC(i) - current rm(1) c parameters C RMD(i) - current rm(2) d parameters C RMCI(i) - inertial contribution from c terms C RMDI(i) - inertial contribution from d terms C C C 'CART' definitions for molecule with NATMS atoms are passed over in C NCART and CARTDF and the calculated principal moments of inertia and C rotational constants are returned in PMI. Data communication with the C calling segment is via the COMMON blocks /MOMI/ and /INPUT/ C C............................................................................ C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=50,maxiso=30,maxcon=3*maxiso,maxrep=20, * maxchg=10,maxato=300, conv=505379.01D0) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /ANCIL/RMC(6),RMD(3),rmci(6),rmdi(3) COMMON /coords/CORDP(MAXATO,3) COMMON /INPUT/CARTRF(MAXATO,4),NPARF(maxpar,3), * NREP(maxpar,maxrep),VPARF(maxpar),NCOBS(MAXCON,2), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2), * VCHANG(maxiso,maxchg),NISOT,NPAR,NOBSC COMMON /HBLK/AA/SBLK/EIGINV REAL*8 CORD(MAXATO,3),VBA(3),TRANS(3,3),PRC(3),RJA(3),VCA(3), 1 A(55),EIGVEC(100),EIGINV(3,3),ABACK(3,3),aa(3,3),aar(3,3) CON=0.174532925199D-1 C C...Initial cartesian coordinates from the CART scheme C CORD(1,1)=0.0D0 CORD(1,2)=0.0D0 CORD(1,3)=0.0D0 N1=1 CORD(2,1)=CARTDF(2,1) CORD(2,2)=0.0D0 CORD(2,3)=0.0D0 IF(CARTDF(3,2))5555,118,5555 118 CS=-0.333333333333D0 SS= 0.942809041582D0 GO TO 121 5555 CS=DCOS(CON*CARTDF(3,2)) SS=DSIN(CON*CARTDF(3,2)) NA=NCART(3,2) 121 IF(NA-N1)128,122,128 122 CORD(3,1)=CORD(NA,1)+CARTDF(3,1)*CS GO TO 129 128 CORD(3,1)=CORD(NA,1)-CARTDF(3,1)*CS 129 CORD(3,2)=CARTDF(3,1)*SS CORD(3,3)=0.0D0 C IF(NATMS.lt.4)goto 5557 DO 5556 I=4,NATMS NO=NCART(I,1) NA=NCART(I,2) NB=NCART(I,3) NC=NCART(I,4) IF(CARTDF(I,2))133,132,133 132 CS=-0.33333333333D0 SS= 0.9428090416D0 GO TO 5558 133 CS=DCOS(CON*CARTDF(I,2)) SS=DSIN(CON*CARTDF(I,2)) 5558 DSQ=0.0D0 DO 138 M=1,3 VBA(M)=CORD(NB,M)-CORD(NA,M) VCA(M)=CORD(NC,M)-CORD(NA,M) DSQ=DSQ+VBA(M)**2 138 CONTINUE RAB=DSQRT(DSQ) SCALE=0.0D0 DO 142 M=1,3 TRANS(M,1)=VBA(M)/RAB SCALE=SCALE+TRANS(M,1)*VCA(M) 142 CONTINUE DSQ=0.0D0 DO 146 M=1,3 RJA(M)=VCA(M)-SCALE*TRANS(M,1) DSQ=DSQ+RJA(M)**2 146 CONTINUE RAJ=DSQRT(DSQ) DO 148 M=1,3 TRANS(M,2)=RJA(M)/RAJ 148 CONTINUE TRANS(1,3)=TRANS(2,1)*TRANS(3,2)-TRANS(3,1)*TRANS(2,2) TRANS(2,3)=TRANS(3,1)*TRANS(1,2)-TRANS(1,1)*TRANS(3,2) TRANS(3,3)=TRANS(1,1)*TRANS(2,2)-TRANS(2,1)*TRANS(1,2) PRC(1)=CARTDF(I,1)*CS PRC(2)=CARTDF(I,1)*SS*DCOS(CON*CARTDF(I,3)) PRC(3)=CARTDF(I,1)*SS*DSIN(CON*CARTDF(I,3)) DO 5556 M=1,3 CORD(NO,M)=CORD(NA,M) DO 5556 K=1,3 CORD(NO,M)=CORD(NO,M)+TRANS(M,K)*PRC(K) 5556 CONTINUE C C...Transform to centre of mass coordinates (using masses of the parent) C 5557 WW=0.0D0 DO 991 I=1,NATMS WW=WW+CARTRF(I,4) 991 CONTINUE FMX=0.0D0 FMY=0.0D0 FMZ=0.0D0 DO 992 N=1,NATMS FMX=FMX+CARTRF(N,4)*CORD(N,1) FMY=FMY+CARTRF(N,4)*CORD(N,2) FMZ=FMZ+CARTRF(N,4)*CORD(N,3) 992 CONTINUE DO 993 N=1,NATMS CORD(N,1)=CORD(N,1)-FMX/WW CORD(N,2)=CORD(N,2)-FMY/WW CORD(N,3)=CORD(N,3)-FMZ/WW 993 CONTINUE C C C Inertial tensor: CART (parent) substituted isotopomer C C 1, = xx = aa C 2,3 xy yy ab bb C 4,5,6 xz yz zz ac bc cc C C symmetric matrix notation: C C corresponding to 1,2,3 for a,b,c - 1,2,3,4,5,6 for a,ab,b,ac,bc,c C after EIGEN - 1,2,3,4,5,6 for c,bc,b,ac,ab,a C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...If this is not the parent isotopomer then: C C 1/ evaluate coordinates of the corresponding parent C 2/ use coordinates of parent to set up I and evaluate its rotation C angles relative to the parent C 3/ determine the I**(1/2) tensor for this isotopomer in inertial C axes of the parent C 4/ add c- and d-elements C 5/ diagonalise the resulting I preserving eigenvalue positions C if(isotop.gt.1)then C C...set up I with masses of the parent and evaluate the parent axes C ix=1 May2009 iy=2 iz=3 418 DO 412 NN=1,6 A(NN)=0.0D0 412 CONTINUE DO 414 N=1,NATMS A(1)=A(1)+CARTRF(N,4)*(CORD(N,2)**2+CORD(N,3)**2) A(2)=A(2)-CARTRF(N,4)* CORD(N,1)* CORD(N,2) A(3)=A(3)+CARTRF(N,4)*(CORD(N,1)**2+CORD(N,3)**2) A(4)=A(4)-CARTRF(N,4)* CORD(N,1)* CORD(N,3) A(5)=A(5)-CARTRF(N,4)* CORD(N,2)* CORD(N,3) A(6)=A(6)+CARTRF(N,4)*(CORD(N,1)**2+CORD(N,2)**2) 414 CONTINUE C CALL EIGEN(A,EIGVEC,3,0) <------- C EIGINV(1,1)=EIGVEC(5)*EIGVEC(3)-EIGVEC(6)*EIGVEC(2) EIGINV(1,2)=EIGVEC(6)*EIGVEC(1)-EIGVEC(4)*EIGVEC(3) EIGINV(1,3)=EIGVEC(4)*EIGVEC(2)-EIGVEC(5)*EIGVEC(1) EIGINV(2,1)=EIGVEC(2)*EIGVEC(9)-EIGVEC(8)*EIGVEC(3) EIGINV(2,2)=EIGVEC(7)*EIGVEC(3)-EIGVEC(9)*EIGVEC(1) EIGINV(2,3)=EIGVEC(8)*EIGVEC(1)-EIGVEC(7)*EIGVEC(2) EIGINV(3,1)=EIGVEC(8)*EIGVEC(6)-EIGVEC(5)*EIGVEC(9) EIGINV(3,2)=EIGVEC(9)*EIGVEC(4)-EIGVEC(7)*EIGVEC(6) EIGINV(3,3)=EIGVEC(7)*EIGVEC(5)-EIGVEC(8)*EIGVEC(4) WW=EIGVEC(7)*EIGINV(1,1)+EIGVEC(4)*EIGINV(2,1) 1 +EIGVEC(1)*EIGINV(3,1) DO 415 N=1,3 DO 415 NN=1,3 EIGINV(N,NN)=EIGINV(N,NN)/WW 415 CONTINUE c DO 416 N=1,NATMS CORDP(N,IX)=CORD(N,1)*EIGINV(1,1)+CORD(N,2)*EIGINV(1,2)+ 1 CORD(N,3)*EIGINV(1,3) CORDP(N,IY)=CORD(N,1)*EIGINV(2,1)+CORD(N,2)*EIGINV(2,2)+ 1 CORD(N,3)*EIGINV(2,3) CORDP(N,IZ)=CORD(N,1)*EIGINV(3,1)+CORD(N,2)*EIGINV(3,2)+ 1 CORD(N,3)*EIGINV(3,3) 416 CONTINUE C C...now set up I in coordinates of the parent that have just been evaluated C (first transferring to centre of mass) C WW=0.0D0 DO 891 I=1,NATMS WW=WW+CARTDF(I,4) 891 CONTINUE FMX=0.0D0 FMY=0.0D0 FMZ=0.0D0 DO 892 N=1,NATMS FMX=FMX+CARTDF(N,4)*CORDP(N,1) FMY=FMY+CARTDF(N,4)*CORDP(N,2) FMZ=FMZ+CARTDF(N,4)*CORDP(N,3) 892 CONTINUE DO 893 N=1,NATMS CORDP(N,1)=CORDP(N,1)-FMX/WW CORDP(N,2)=CORDP(N,2)-FMY/WW CORDP(N,3)=CORDP(N,3)-FMZ/WW 893 CONTINUE C if(ncall.eq.0.and.idump.eq.1)then write(3,'(1x/1x,39(2h -)/30x, * ''Coordinates in parent axes for isotope'', * I3,'':'')')isotop do 515 n=1,natms write(3,'(30x,3f14.8,f14.8)') * (cordp(n,nn),nn=1,3),cartdf(n,4) 515 continue endif C DO 994 N=1,3 do 994 nn=1,3 AA(N,NN)=0.0D0 994 CONTINUE C DO 995 N=1,NATMS AA(1,1)=AA(1,1)+CARTDF(N,4)*(CORDP(N,2)**2+CORDP(N,3)**2) AA(1,2)=AA(1,2)-CARTDF(N,4)* CORDP(N,1)* CORDP(N,2) AA(2,2)=AA(2,2)+CARTDF(N,4)*(CORDP(N,1)**2+CORDP(N,3)**2) AA(1,3)=AA(1,3)-CARTDF(N,4)* CORDP(N,1)* CORDP(N,3) AA(2,3)=AA(2,3)-CARTDF(N,4)* CORDP(N,2)* CORDP(N,3) AA(3,3)=AA(3,3)+CARTDF(N,4)*(CORDP(N,1)**2+CORDP(N,2)**2) 995 CONTINUE AA(2,1)=AA(1,2) AA(3,1)=AA(1,3) AA(3,2)=AA(2,3) C DO 996 N=1,3 DO 996 NN=1,3 ABACK(N,NN)=AA(N,NN) 996 CONTINUE C if(ncall.eq.0.and.idump.eq.1)then write(3,'( * 30x,''I in parent axes for isotope'',i3,'':''/30x, * 3f14.8/30x,3f14.8/30x,3f14.8)')isotop, * ((ABACK(n,NN),nn=1,3),n=1,3) endif C C...diagonalise, preserving eigenvalue positions. Resulting eigenvectors C describe rotation from axes of the parent C CALL HDIAG(3,0,0) <------- C C...Safeguard against oblate axis switching C if(AA(1,1).gt.AA(2,2)+0.01d0)then May2009 if(ix.eq.1.and.iy.eq.2)then ix=2 iy=1 goto 418 endif endif C C...Safeguard against prolate axis switching C if(AA(2,2).gt.AA(3,3)+0.01d0)then May2009 if(iy.eq.2.and.iz.eq.3)then iy=3 iz=2 goto 418 endif endif C if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''Eigenvectors of I for isotope'',i3)')isotop do 500 n=1,3 write(3,'(30x,3f14.8)')(eiginv(n,nn),nn=1,3) 500 continue write(3,'(30x,''Eigenvalues of I:''/30x,3f14.8)') * AA(1,1),aa(2,2),aa(3,3) endif C C...I**(1/2) matrix of the isotope in axes of the parent C do 1000 n=1,3 do 1000 nn=1,3 if(n.ne.nn)aa(n,nn)=0.d0 aar(n,nn)=0.d0 1000 continue AA(1,1)=dsqrt(AA(1,1)) AA(2,2)=dsqrt(AA(2,2)) AA(3,3)=dsqrt(AA(3,3)) do 1001 n=1,3 do 1001 nn=1,3 do 1002 nnn=1,3 aar(n,nn)=aar(n,nn)+eiginv(n,nnn)*eiginv(nn,nnn) * *aa(nnn,nnn) 1002 continue 1001 continue C if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''I**(1/2) matrix for isotope'',i3)')isotop do 1500 n=1,3 write(3,'(30x,3f14.8)')(aar(n,nn),nn=1,3) 1500 continue endif C DO 999 N=1,3 DO 999 NN=1,3 AA(N,NN)=ABACK(N,NN) 999 CONTINUE C C...rm() c contribution C RMCI(1)=AAR(1,1)*RMC(1) RMCI(2)=AAR(2,2)*RMC(2) RMCI(3)=AAR(3,3)*RMC(3) RMCI(4)=AAR(1,2)*RMC(4) RMCI(5)=AAR(1,3)*RMC(5) RMCI(6)=AAR(2,3)*RMC(6) AA(1,1)=AA(1,1)+RMCI(1) AA(2,2)=AA(2,2)+RMCI(2) AA(3,3)=AA(3,3)+RMCI(3) AA(1,2)=AA(1,2)+RMCI(4) AA(1,3)=AA(1,3)+RMCI(5) AA(2,3)=AA(2,3)+RMCI(6) AA(2,1)=AA(1,2) AA(3,1)=AA(1,3) AA(3,2)=AA(2,3) C C...rm() d contribution C www=1.d0 swww=0.d0 DO 312 N=1,NATMS atmass=cartdf(n,4) if(atmass.gt.0.5d0)then www=www*ATMASS swww=swww+ATMASS endif 312 CONTINUE www=(www/swww)**(1.d0/real(2*(NATMS-NDATMS)-2)) RMDI(1)=www*RMD(1) RMDI(2)=www*RMD(2) RMDI(3)=www*RMD(3) AA(1,1)=AA(1,1)+RMDI(1) AA(2,2)=AA(2,2)+RMDI(2) AA(3,3)=AA(3,3)+RMDI(3) C if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''I+c+d in parent axes for isotope'',i3/30x, * 3f14.8/30x,3f14.8/30x,3f14.8)')isotop, * ((AA(n,NN),nn=1,3),n=1,3) endif C CALL HDIAG(3,0,0) <------- RIA=AA(1,1) RIB=AA(2,2) RIC=AA(3,3) C endif C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...If this is the parent isotopomer then: C C 1/ Use initial cm cartesians to set up and diagonalise I C 2/ Evaluate principal coordinates of the parent for later use by the C rm(2) rotation of axes C 3/ Add c- and d- terms (in this case only the diagonal c-terms apply) C if(isotop.eq.1)then if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''Pre-diagonalisation coords of parent:'')') do 710 n=1,natms write(3,'(30x,3f14.8,f14.8)') * (cord(n,nn),nn=1,3),cartdf(n,4) 710 continue endif C DO 1994 N=1,6 A(N)=0.0D0 1994 CONTINUE DO 1995 N=1,NATMS A(1)=A(1)+CARTDF(N,4)*(CORD(N,2)**2+CORD(N,3)**2) A(2)=A(2)-CARTDF(N,4)* CORD(N,1)* CORD(N,2) A(3)=A(3)+CARTDF(N,4)*(CORD(N,1)**2+CORD(N,3)**2) A(4)=A(4)-CARTDF(N,4)* CORD(N,1)* CORD(N,3) A(5)=A(5)-CARTDF(N,4)* CORD(N,2)* CORD(N,3) A(6)=A(6)+CARTDF(N,4)*(CORD(N,1)**2+CORD(N,2)**2) 1995 CONTINUE C if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''I for the parent:''/30x, * f14.8/30x,2f14.8/30x,3f14.8)')(A(n),n=1,6) endif C CALL EIGEN(A,EIGVEC,3,0) <------- c if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''Eigenvalues of I for the parent:'' * /30x,3f14.8)')A(1),a(3),a(6) endif c RIA=A(6) RIB=A(3) RIC=A(1) C EIGINV(1,1)=EIGVEC(5)*EIGVEC(3)-EIGVEC(6)*EIGVEC(2) EIGINV(1,2)=EIGVEC(6)*EIGVEC(1)-EIGVEC(4)*EIGVEC(3) EIGINV(1,3)=EIGVEC(4)*EIGVEC(2)-EIGVEC(5)*EIGVEC(1) EIGINV(2,1)=EIGVEC(2)*EIGVEC(9)-EIGVEC(8)*EIGVEC(3) EIGINV(2,2)=EIGVEC(7)*EIGVEC(3)-EIGVEC(9)*EIGVEC(1) EIGINV(2,3)=EIGVEC(8)*EIGVEC(1)-EIGVEC(7)*EIGVEC(2) EIGINV(3,1)=EIGVEC(8)*EIGVEC(6)-EIGVEC(5)*EIGVEC(9) EIGINV(3,2)=EIGVEC(9)*EIGVEC(4)-EIGVEC(7)*EIGVEC(6) EIGINV(3,3)=EIGVEC(7)*EIGVEC(5)-EIGVEC(8)*EIGVEC(4) WW=EIGVEC(7)*EIGINV(1,1)+EIGVEC(4)*EIGINV(2,1) 1 +EIGVEC(1)*EIGINV(3,1) DO 3060 N=1,3 DO 3060 NN=1,3 EIGINV(N,NN)=EIGINV(N,NN)/WW 3060 CONTINUE DO 998 N=1,NATMS CORDP(N,1)=CORD(N,1)*EIGINV(1,1)+CORD(N,2)*EIGINV(1,2)+ 1 CORD(N,3)*EIGINV(1,3) CORDP(N,2)=CORD(N,1)*EIGINV(2,1)+CORD(N,2)*EIGINV(2,2)+ 1 CORD(N,3)*EIGINV(2,3) CORDP(N,3)=CORD(N,1)*EIGINV(3,1)+CORD(N,2)*EIGINV(3,2)+ 1 CORD(N,3)*EIGINV(3,3) 998 CONTINUE C if(ncall.eq.0.and.idump.eq.1)then write(3,'(30x,''Coordinates of parent:'')') do 510 n=1,natms write(3,'(30x,3f14.8,f14.8)') * (cordp(n,nn),nn=1,3),cartdf(n,4) 510 continue endif C C...rm() c contribution (only diagonal terms in this case) C RMCI(1)=DSQRT(RIA)*RMC(1) RMCI(2)=DSQRT(RIB)*RMC(2) RMCI(3)=DSQRT(RIC)*RMC(3) RMCI(4)=0.D0 RMCI(5)=0.D0 RMCI(6)=0.D0 RIA=RIA+RMCI(1) RIB=RIB+RMCI(2) RIC=RIC+RMCI(3) C C...rm() d contribution C www=1.d0 swww=0.d0 DO 301 N=1,NATMS atmass=cartdf(n,4) if(atmass.gt.0.5d0)then www=www*ATMASS swww=swww+ATMASS endif 301 CONTINUE www=(www/swww)**(1.d0/real(2*(NATMS-NDATMS)-2)) RMDI(1)=www*RMD(1) RMDI(2)=www*RMD(2) RMDI(3)=www*RMD(3) RIA=RIA+RMDI(1) RIB=RIB+RMDI(2) RIC=RIC+RMDI(3) endif C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...Calculated moments of inertia and planar moments C PMI(1,1)=RIA PMI(2,1)=RIB PMI(3,1)=RIC C PMI(4,1)=4.0d0*RIB*RIC/(RIB+RIC) PMI(5,1)=4.0d0*RIA*RIB/(RIA+RIB) PMI(6,1)=0.5d0*(RIB+RIC-RIA) PMI(7,1)=0.5d0*(RIA+RIC-RIB) PMI(8,1)=0.5d0*(RIA+RIB-RIC) C C...Calculated rotational constants C PMI(1,2)=conv/PMI(1,1) PMI(2,2)=conv/PMI(2,1) PMI(3,2)=conv/PMI(3,1) C PMI(4,2)=conv*4.d0/PMI(4,1) PMI(5,2)=conv*4.d0/PMI(5,1) if(PMI(6,1).ne.0.d0)then PMI(6,2)=conv/PMI(6,1) else PMI(6,2)=1.D+15 endif if(PMI(7,1).ne.0.d0)then PMI(7,2)=conv/PMI(7,1) else PMI(7,2)=1.D+15 endif if(PMI(8,1).ne.0.d0)then PMI(8,2)=conv/PMI(8,1) else PMI(8,2)=1.D+15 endif C RETURN END C C_____________________________________________________________________________ C SUBROUTINE EIGEN(A,R,N,MV) C C MATRIX DIAGONALIZATION ROUTINE FOR SYMMETRIC MATRICES FROM THE C 'SSP' LIBRARY C C - Jacobi type algorithm C - matrix to be diagonalized has to be in symmetric storage mode C (A and R are presently dimensioned for 10x10 matrices) C DIMENSION A(55),R(100) DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX, 1 COSX2,SINCS,RANGE,DSQRT,DABS C C C GENERATE IDENTITY MATRIX C 5 RANGE=1.0D-12 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0.0d0 IF(I-J) 20,15,20 15 R(IJ)=1.0d0 20 CONTINUE C C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) C 25 ANORM=0.0d0 DO 35 I=1,N DO 35 J=I,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM) 165,165,40 40 ANORM=1.41421356237d0*DSQRT(ANORM) ANRMX=ANORM*RANGE/FLOAT(N) C C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR C IND=0 THR=ANORM 45 THR=THR/FLOAT(N) 50 L=1 55 M=L+1 C C COMPUTE SIN AND COS C 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(DABS(A(LM))-THR) 130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=0.5d0*(A(LL)-A(MM)) 68 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/DSQRT(2.0d0*(1.0d0+(DSQRT(1.0d0-Y*Y)))) SINX2=SINX*SINX 78 COSX=DSQRT(1.0d0-SINX2) COSX2=COSX*COSX SINCS =SINX*COSX C C ROTATE L AND M COLUMNS C ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GO TO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GO TO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.0d0*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X C C TESTS FOR COMPLETION C C TEST FOR M = LAST COLUMN C 130 IF(M-N) 135,140,135 135 M=M+1 GO TO 60 C C TEST FOR L = SECOND FROM LAST COLUMN C 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GO TO 55 150 IF(IND-1) 160,155,160 155 IND=0 GO TO 50 C C COMPARE THRESHOLD WITH FINAL NORM C 160 IF(THR-ANRMX) 165,165,45 C C SORT EIGENVALUES AND EIGENVECTORS C 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HDIAG(N,IEGEN,IORD) C C A=MATRIX TO BE DIAGONALIZED C N=DIMENSION OF SUBMATRIX TO BE DIAGONALIZED C IEGEN=0 EIGENVALUES AND EIGENVECTORS C =1 EIGENVALUES ONLY C IORD=0 NO ORDERING OF EIGENVALUES OR EIGENVECTORS (ORDER IN = ORDER OUT) C =1 ORDERING BY SIZE OF EIGENVALUES C EIVR=TRANSFORMATION MATRIX (MATRIX OF EIGENVECTORS) C C THIS ROUTINE USES A VARIABLE THRESHOLD JACOBI METHOD C IMPLICIT REAL*8(A-H,O-Z) PARAMETER (MDIM=3) COMMON /HBLK/A/SBLK/EIVR DIMENSION A(MDIM,MDIM),EIVR(MDIM,MDIM) c c write(*,'('' starting '',i4,''x'',i4)')n,n IF(N.GT.0)GOTO 1 EIVR(1,1)=1.0 RETURN 1 IF(IEGEN.GT.0)GOTO 102 DO 101 J=1,N DO 100 I=1,N EIVR(I,J)=0.0d0 100 CONTINUE EIVR(J,J)=1.0d0 101 CONTINUE C C FIND THE ABSOLUTELY LARGEST ELEMENT OF A C 102 ATOP=0.0D0 DO 111 I=1,N DO 111 J=1,N IF(ATOP.GE.DABS(A(I,J)))GOTO 111 ATOP=DABS(A(I,J)) 111 CONTINUE IF(ATOP)109,109,113 109 RETURN C C CALCULATE THE STOPPING CRITERION -- DSTOP C 113 AVGF=(N*(N-1))*0.55D0 D=0.0D0 DO 114 JJ=2,N DO 114 II=2,JJ S=A(II-1,JJ)/ATOP D=S*S+D 114 CONTINUE DSTOP=(1.0D-06)*D C C CALCULATE THE THRESHOLD -- THRSH C THRSH=DSQRT(D/AVGF)*ATOP C C START A SWEEP C 115 IFLAG=0 DO 130 JCOL=2,N JCOL1=JCOL-1 DO 130 IROW=1,JCOL1 AIJ=A(IROW,JCOL) C C COMPARE THE OFF DIAGONAL ELEMENT WITH THRSH C IF(DABS(AIJ).LE.THRSH)GOTO 130 AII=A(IROW,IROW) AJJ=A(JCOL,JCOL) S=AJJ-AII C C CHECK TO SEE IF THE CHOSEN ROTATION IS LESS THAN THE ROUNDING ERROR C IF SO, THEN DO NOT ROTATE C IF(DABS(AIJ).LE.(1.0D-09*DABS(S)))GOTO 130 IFLAG=1 C C IF THE ROTATION IS VERY CLOSE TO 45 DEGREES, SET SIN AND COS C TO 1/(ROOT 2) C IF((1.0D-10*DABS(AIJ)).LT.DABS(S))GOTO 116 S=0.70710678118654752D0 C=S GOTO 120 C C CALCULATE SIN AND COS FOR ROTATION THAT IS NOT VERY CLOSE C TO 45 DEGREES C 116 T=AIJ/S S=0.25D0/DSQRT(0.25D0+T*T) C C COS=C, SIN=S C C=DSQRT(0.5D0+S) S=2.0D0*T*S/C C C CALCULATION OF THE NEW ELEMENTS OF MATRIX A C 120 DO 121 I=1,IROW T=A(I,IROW) U=A(I,JCOL) A(I,IROW)=C*T-S*U A(I,JCOL)=S*T+C*U 121 CONTINUE I2=IROW+2 IF(I2.GT.JCOL)GOTO 123 CONTINUE DO 122 I=I2,JCOL T=A(I-1,JCOL) U=A(IROW,I-1) A(I-1,JCOL)=S*U+C*T A(IROW,I-1)=C*U-S*T 122 CONTINUE 123 A(JCOL,JCOL)=S*AIJ+C*AJJ A(IROW,IROW)=C*A(IROW,IROW)-S*(C*AIJ-S*AJJ) DO 124 J=JCOL,N T=A(IROW,J) U=A(JCOL,J) A(IROW,J)=C*T-S*U 124 A(JCOL,J)=S*T+C*U C C ROTATION COMPLETED C SEE IF EIGENVECTORS REQUIRED C IF(IEGEN.GT.0)GOTO 126 DO 125 I=1,N T=EIVR(I,IROW) EIVR(I,IROW)=SNGL(C*T-EIVR(I,JCOL)*S) EIVR(I,JCOL)=SNGL(S*T+EIVR(I,JCOL)*C) 125 CONTINUE C C CALCULATE THE NEW NORM D AND COMPARE WITH DSTOP C 126 CONTINUE S=AIJ/ATOP D=D-S*S IF(D.GE.DSTOP)GOTO 129 C C RECALCULATE DSTOP AND THRSH TO DISCARD ROUNDING ERRORS C D=0.0D0 DO 128 JJ=2,N DO 128 II=2,JJ S=A(II-1,JJ)/ATOP D=S*S+D 128 CONTINUE DSTOP=(1.0D-06)*D 129 THRSH=DSQRT(D/AVGF)*ATOP 130 CONTINUE IF(IFLAG.NE.0)GOTO 115 C C ARRANGE THE EIGENVALUES IN THE ORDER OF INCREASING ENERGY C ARRANGE THE EIGENVECTORS IN THE SAME ORDER C IF(IORD.EQ.0)RETURN NU=N DO 11 I=1,N IF(I.GE.NU)RETURN AMIN=A(I,I) DO 10 J=I,NU IF(A(J,J).GE.AMIN)GOTO 10 C C IF IEGEN IS -1 EXCLUDE UNCONVERGED EIGENVALUES FROM THIS ORDERING C TE=ABS(EIVR(N,J))+ABS(EIVR(N-1,J)) IF((TE.GT.0.05D0).AND.(IEGEN.EQ.-1))GOTO 15 II=I AMIN=A(J,J) A(J,J)=A(I,I) A(I,I)=AMIN 16 DO 12 K=1,N TEMP=EIVR(K,II) EIVR(K,II)=EIVR(K,J) EIVR(K,J)=SNGL(TEMP) 12 CONTINUE GOTO 10 15 AM=A(J,J) A(J,J)=A(NU,NU) A(NU,NU)=AM II=NU NU=NU-1 GOTO 16 10 CONTINUE 11 CONTINUE C RETURN END C C_____________________________________________________________________________ C C SUBROUTINE MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT, 1 COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA) C C............................................................................. C C Levenberg-Marquardt method, attempting to reduce the value CHISQ of a fit C between a set of NDATA points X(I),Y(I) with individual standard devia- C tions SIG(I), and a nonlinear function dependent on MA coefficients A. C C C On the first call of MRQMIN provide an initial guess for the parameters C A, and set ALAMDA<0 for initialization (which then sets ALAMDA=0.001). C If a step succeeds CHISQ becomes smaller and ALAMDA decreases C by a factor of 10. C If a step fails ALAMDA grows by a factor of 10. You must call this C routine repeatedly until convergence is achieved. Then make one final C call with ALAMDA=0., so that COVAR(I,J) returns the covariance matrix, C and ALPHA(I,J) the curvature matrix. C C NOTE: This is double precision adaptation of the routine set from Num.Rec. C and there is also an extra CALL SETUP statement in MRQCOF to C work with STRFIT - these and other modifications identified C with ...STRMOD C............................................................................. C C ROUTINES REQUIRED: C C MRQCOF, COVSRT, GAUSSJ - called by MRQMIN C FUNCS(user supplied) - called by MRQCOF C C C X,Y,SIG - vectors of length NDATA containin X,Y pairs and deviations C on the Y values C A - vector of length MA containing the constants to be used for C evaluations of Y.calc (MFIT of which are to be fitted) C LISTA - vector of length MFIT carrying pointers to the constants in C A which are to be fitted C COVAR,ALPHA - square arrays with physical dimension NCA (which is to be C greater or equal to MFIT) used as workspace. Routine C COVSRT will return the covariance matrix in COVAR. C CHISQ - minimisation parameter chi-square C FUNCS - user supplied subroutine FUNCS(X,A,YFIT,DYDA,MA) that C evaluates the calculated value of Y: YFIT at X, and the C derivatives of Y with respect to the fitting parameters A C at X. C ALAMDA - parameter to monitor the prograss of fit, to be set to any C -ve value for initialisation. This will then be set to 0.001 C and on following cycles decreased by a factor of 10 is chi- C square decreases and increased by a factor of 10 if CHISQ C increases. ALAMDA=0 returns the covariance matrix COVAR C and the curvature (Hessian) matrix ALPHA C............................................................................. C C...MMAX defines the largest number of parameters C IMPLICIT REAL*8 (A-H,O-Z) ...STRMOD EXTERNAL FUNCS ...STRMOD PARAMETER (maxpar=50, MMAX=maxpar) ...STRMOD COMMON /COV/COVTMP(maxpar,maxpar) ...STRMOD DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MFIT), 1 COVAR(NCA,NCA),ALPHA(NCA,NCA),ATRY(MMAX),BETA(MMAX),DA(MMAX) c c...initialization c IF(ALAMDA.LT.0.d0)THEN KK=MFIT+1 c c...does LISTA contain a proper permutation of the coefficients? DO 12 J=1,MA IHIT=0 DO 11 K=1,MFIT IF(LISTA(K).EQ.J)IHIT=IHIT+1 11 CONTINUE IF (IHIT.EQ.0) THEN LISTA(KK)=J KK=KK+1 ELSE IF (IHIT.GT.1) THEN PAUSE 'Improper permutation in LISTA' ENDIF 12 CONTINUE IF (KK.NE.(MA+1)) PAUSE 'Improper permutation in LISTA' ALAMDA=0.001d0 CALL MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NCA,CHISQ, 1 FUNCS) <------- OCHISQ=CHISQ DO 13 J=1,MA ATRY(J)=A(J) 13 CONTINUE ENDIF c c...alter linearized fitting matrix, by augmenting diagonal elements c DO 15 J=1,MFIT DO 14 K=1,MFIT COVAR(J,K)=ALPHA(J,K) 14 CONTINUE COVAR(J,J)=ALPHA(J,J)*(1.+ALAMDA) DA(J)=BETA(J) 15 CONTINUE c c...matrix solution c CALL GAUSSJ(COVAR,MFIT,NCA,DA,1,1) <------- c c...once converged evaluate covariance matrix with ALAMDA=0 c (but if ALAMDA not zero then transfer covariance matrix to ...STRMOD c array COVTMP for passing over to calling program for evaluation ...STRMOD c of errors for each cycle) ...STRMOD c IF(ALAMDA.NE.0.d0)THEN ...STRMOD DO 20 IZA=1,MA ...STRMOD DO 20 IZB=1,MA ...STRMOD COVTMP(IZA,IZB)=COVAR(IZA,IZB) ...STRMOD 20 CONTINUE ...STRMOD ENDIF IF(ALAMDA.EQ.0.d0)THEN CALL COVSRT(COVAR,NCA,MA,LISTA,MFIT) <------- RETURN ENDIF c c...did the trial succeed? c DO 16 J=1,MFIT ATRY(LISTA(J))=ATRY(LISTA(J))+DA(J) 16 CONTINUE CALL MRQCOF(X,Y,SIG,NDATA,ATRY,MA,LISTA,MFIT,COVAR,DA,NCA,CHISQ, 1 FUNCS) <------- c c...success, accept the new solution c IF(CHISQ.LT.OCHISQ)THEN ALAMDA=0.1d0*ALAMDA OCHISQ=CHISQ DO 18 J=1,MFIT DO 17 K=1,MFIT ALPHA(J,K)=COVAR(J,K) 17 CONTINUE BETA(J)=DA(J) A(LISTA(J))=ATRY(LISTA(J)) 18 CONTINUE c c...failure, increase ALAMDA and return c ELSE ALAMDA=10.d0*ALAMDA C C the following has been commented out to avoid termination of iterations C CHISQ=OCHISQ ...STRMOD ENDIF c RETURN END c c_____________________________________________________________________________ c SUBROUTINE MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NALP, 1 CHISQ,FUNCS) C C Routine used by MRQMIN to evaluate the linearized fitting matrix ALPHA, C and vector BETA (equation 14.4.8 Num.Rec.) C C IMPLICIT REAL*8 (A-H,O-Z) ...STRMOD PARAMETER (maxpar=50, MMAX=maxpar) ...STRMOD DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),ALPHA(NALP,NALP),BETA(MA), * DYDA(MMAX),LISTA(MFIT),A(MA) c c...initialize (symmetric) ALPHA,BETA c DO 12 J=1,MFIT DO 11 K=1,J ALPHA(J,K)=0.d0 11 CONTINUE BETA(J)=0.d0 12 CONTINUE CHISQ=0. c c...summation loop over all data (the call to SETUP is the only modification c of the Numerical Recipes version for STRFIT) c CALL SETUP(A,MA) <--...STRMOD C DO 15 I=1,NDATA CALL FUNCS(X(I),A,YMOD,DYDA,MA) <--FUNCS SIG2I=1.d0/(SIG(I)*SIG(I)) DY=Y(I)-YMOD DO 14 J=1,MFIT WT=DYDA(LISTA(J))*SIG2I DO 13 K=1,J ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K)) 13 CONTINUE BETA(J)=BETA(J)+DY*WT 14 CONTINUE c c...find CHISQ CHISQ=CHISQ+DY*DY*SIG2I 15 CONTINUE c c...fill in the symmetric side c DO 17 J=2,MFIT DO 16 K=1,J-1 ALPHA(K,J)=ALPHA(J,K) 16 CONTINUE 17 CONTINUE c RETURN END c c_____________________________________________________________________________ c SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) C C Routine for solving a set of linear equation by the Gauss-Jordan elimina- C tion with full pivoting. C C A is an input matrix of N by N elements, stored in an array of physical C dimensions NP by NP. C B is an input matrix of N by M containing the M right-hand side vectors, C stored in an array of physical dimensions NP by MP. C On output, A is replaced by its matrix inverse, and B is replaced by the C corresponding set of solution vectors C C IMPLICIT REAL*8 (A-H,O-Z) ...STRMOD PARAMETER (NMAX=50) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) c DO 11 J=1,N IPIV(J)=0 11 CONTINUE c c...this is the main loop over the columns to be reduced c DO 22 I=1,N BIG=0.d0 c c...this is the outer loop of the search for the pivot element c DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'Singular matrix' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 c c...we now have the pivot element, so we interchange rows, if needed to put c the pivot on the diagonal. The columns are not physically interchanged, c only relabeled: c IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.d0) PAUSE 'Singular matrix.' PIVINV=1.d0/A(ICOL,ICOL) A(ICOL,ICOL)=1.d0 DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0.d0 DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE c c...this is the end of the main loop over columns of the reduction. c It only remains to unscramble the solution in view of the column c interchanges. c DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE c RETURN END c c_____________________________________________________________________________ c SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT) C C Given the covariance matrix COVAR of a fit for MFIT of MA total parameters C and their ordering LISTA(I), repack the covariance matrix to the true C order of the parameters. C Elements associated with fixed parameters will be zero. NCVM is the C physical dimension of COVAR C C IMPLICIT REAL*8 (A-H,O-Z) ...STRMOD DIMENSION COVAR(NCVM,NCVM),LISTA(MFIT) c c...zero all elements below diagonal c DO 12 J=1,MA-1 DO 11 I=J+1,MA COVAR(I,J)=0.d0 11 CONTINUE 12 CONTINUE c c...repack off-diagonal elements of fit into correct locations below diagonal c DO 14 I=1,MFIT-1 DO 13 J=I+1,MFIT IF(LISTA(J).GT.LISTA(I)) THEN COVAR(LISTA(J),LISTA(I))=COVAR(I,J) ELSE COVAR(LISTA(I),LISTA(J))=COVAR(I,J) ENDIF 13 CONTINUE 14 CONTINUE c c...temporarily store original diagonal elements in top row, and zero the c diagonal c SWAP=COVAR(1,1) DO 15 J=1,MA COVAR(1,J)=COVAR(J,J) COVAR(J,J)=0.d0 15 CONTINUE COVAR(LISTA(1),LISTA(1))=SWAP c c...now sort elements into proper order on diagonal c DO 16 J=2,MFIT COVAR(LISTA(J),LISTA(J))=COVAR(1,J) 16 CONTINUE c c...finally fill in above diagonal by symmetry c DO 18 J=2,MA DO 17 I=1,J-1 COVAR(I,J)=COVAR(J,I) 17 CONTINUE 18 CONTINUE c RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________ C C C PROGRAM FLOW IN M/P: C | C |__DATAIN read data C | C |__MRQMIN initialisation with ALAMDA<0 C | C |__MRQMIN <--- minimisation runs, first with C | |_______| ALAMDA=0.001 C | loop back if CHISQ not converged C | C |__MRQMIN final call with ALAMDA=0. C | C |__write output files C C C PROGRAM FLOW IN MRQMIN: C | C |__MRQCOF C | |____SETUP zeroth call for moments of inertia C | | |___SETUP1 for the current data C | | | |___MOMIN C | | | |___EIGEN C | | | C | | |___SETUP1 further calls (equal to the number C | | of parameters) for derivatives C | | C | |____FUNCS = LMUSER load up Y and DYDA evaluated above C | | C | |____GAUSSJ solve linear equations C | | C | |____COVSRT tidy up covariance matrix C | C |__MRQCOF second call to establish change of C CHISQ C C_____________________________________________________________________________ C_____________________________________________________________________________