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 it can also C be used to fit rz, rm^(1L), rm^(2L), re^(SE) and other geometries. C C Citation and additional 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 Sciences, 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 basic least-squares engine for fitting r0 internal parameters to C ground state moments of inertia can be used to evaluate also: C C 1/ r.z, re^SE or any other geometry fitted to moments of inertia C modified by user supplied corrective contributions. For example, C in the re^SE geometry (= semi-experimental equilibrium) the Be-B0 C corrections from a CFOUR anharmonic force field calculation can C be used directly. The much smaller rotational g-tensor contributions C can also be used. 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 Reviews dealing (among other topics) with re^SE evaluations: C J.Demaison, Molec.Phys. 105,3109-3138(2007) C C.Puzzarini et al., Int.Rev.Pjys.Chem. 29,273-367(2010) C C Demonstration of the usefulness of rotational g-tensor corrections C in r_e^SE geometry evaluations: C A.G.Csaszar,J.Demaison,H.D.Rudolph, "Equilibrium structures of three- C four-, five-, six-, and seven-membered unsaturated N-containing C heterocycles", J.Phys.Chem.A 119,1731-1746(2015) 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 C - each structural parameter can have up to twenty repetitions in the C molecule C C - up to MAXCHG structural parameters can be different between the parent C and the substituted species C C - non-linear Levenberg-Marquardt least squares fitting algorithm is used C (Numerical Recipes Chapt.14) C 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 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 C - the data set for STRFIT is compatible with PMIFST so the latter can C be used to check the 'CART' definitions C C - The user only needs to specify the generic name MOLNAM on startup, C and the input/output files are then: C C MOLNAM.STF = the input file C C MOLNAM.OUT = the main output file with complete results of the fit C MOLNAM.PMI = the abbreviated output file containing at the top C the connectivity declaration of the fitted structure C for viewing with PMIFST C MOLNAM.XYZ = the fitted geometry in Cartesians for viewing with packages C such as Chimera, JMol, MacMolPlot... C C C ver. 14.VI.2021 ----- 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 15.12.10: optional user parameter descriptors and updated input description C 17.02.12: uncertainties in principal coordinates C 12.12.12: increased dimensioning C 6.03.13: allowance for exclusion from fit of some of the declared constants C 6.05.15: parameter repetitions with changes c 22.05.15: generic i/o file names and addition of .xyz output c 8.10.16: improved diagnostics of errors in the input data file C 25.07.19: electronic correction in r_e^SE evaluation + misc. updates C 8.04.20: some output and documentation changes C 19.05.20: removing code generating "Fortran2018 deleted feature" + other warnings C 14.06.21: clarifying printed deviations and output mods C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C STRFIT 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 C The various lines are: C C line 1: descriptive comment (truncated to 78 characters if longer) c C------------------------------------------------------------------------------- C lines 2-13 - block containing the CART definition of the molecular geometry C------------------------------------------------------------------------------- C C line 2: The number of atoms in the molecule, read in with format=(I5) C Optionally a second number (also in an I5 field) can be included, C and if equal to 1 specifies output of debugging information. C C line 3+: CART DEFINITION LINE: C N,NA,NB,NC,R,A,D,W C C this is read in free format, but N,NA,NB,NC are to be integers, and C R,A,D,W are to be floating point values. There should be as many C such lines as the number of atoms declared in Line 2. C C The structure of the CART line is defined in a separate set of C comments further below. C CART lines declare either real or dummy atoms, and the use of dummy C atoms is encouraged since it may considerably simplify the declaration C of the molecule. 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 C C------------------------------------------------------------------------------- C lines14-16 - block defining the parameters of fit C------------------------------------------------------------------------------- C C line 14: The total number of declared parameters. These parameters can either C be parameters of fit or are to be held fixed at declared 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 bond length C change (mostly for use with deuteration) C C Parameter of fit declaration lines 15+ can be in two different forms, 15a or 15b: C C C line 15a: DEFINITION OF AN INTERNAL COORDINATE: C N,I,J,V,AD C C this is read in with input format=(31x,3i2,f14.0,20x,A20) so that C the numbers have to be EXACTLY in the columns indicated below: C C column: 32 38 51 72 91 C | | | | | C --++--++++++++++++++ cccccccccccccccccccc C Examples: C C atom no., parameter no. 7 1 0 r(C-C) C atom no., parameter no. 7 2 2 CCH (in CH3) C 8 9 C atom no., parameter no. FIX 8 2 0 15.346 fixed at lit.value C C N - the number of atom containing the required internal coordinate in C its CART definition C I - the CART parameter for that atom corresponding to the required internal C (1=bond length, 2=bond angle, 3=dihedral angle) C C J - the optional repeat count declaring whether this parameter C has any additional occurrences: C J = 0 means that the parameter is unique C J > 0 makes it possible to fit 'degenerate' parameters, C occurring for more than one atom (eg. CX bond in a CX3 pyramidal C segment) and the value of J defines the number C of additional occurrences of the parameter (e.g. J=2 for bond C lengths in a methyl group). An additional line C in format =(a31,20i2) should then be inserted immediately C following the current definition line and containing C the numbers of atoms also carrying this parameter. C Up to twenty 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 atoms have to be smaller C than 10 (single digit) C J < 0 makes it possible to fit 'degenerate' parameters with changes. C This feature is designed primarily for dihedral angles in C order to allow implementation of local rotational symmetry based C on a non-zero rotation angle for atom N. C This line is to be followed by |J| lines, each containing C two numbers: integer number of the atom carrying the repeated value, C and the difference of this value relative to that for atom N. C The entries are to start from column 32, but are free-format C otherwise. Negative atom numbers are allowed, and these will C specify a reversed sign value of V (see below). With this method C there is no limitation on the atom number that is present for J>0. C C V - the optional initial value for the declared parameter, C which will override the value in the preceding CART deck C AD - the optional alphanumeric descriptor of this parameter C (up to 20 characters long), which will be echoed in the output. C This has to start in column 72. C FIX - the command to fix this parameter in the fit, details at the end C of 15b below C C C line 15b: DEFINITION OF AN ANCILLARY PARAMETER OF FIT: C N,M,I,V C C this is read in with format=(31x,3i2,f14.0) so that the numbers have C to be EXACTLY in the columns indicated below: C C column: 32 38 51 C | | | C --++--++++++++++++++ C C Examples: C C rm(1) parameter c_a........-4 1 0 0.001 C rm(2) parameter_d_b........-5 2 0 C rm(2) parameter_d_b.....FIX-5 2 0 0.001 C C N - the code (always negative) defining the type of the ancillary C parameter, and also the meaning of M,I,V C C N = -4 Watson's rm(1) c parameter, in which case C M = 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 I = 0 C V = optional initial value for this parameter C C N = -5 Watson's rm(2) d parameter, in which case C M = 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 I = 0 C V = optional initial value for this parameter, which will override C the value in the CART definitions C C N = -6 Laurie correction term for an XH bond, in which case C M = the atom for which this bond is defined C I = the number of additional bonds subject to the same correction: if C this is greater than 0 then an additional C line, format=(a31,20i2) should be inserted immediately C following, with numbers of atoms for which these internals C are defined C V = optional initial value for this parameter C C N = -7 experimental Laurie term defined in terms of a projection angle C M = the atom for which this angle is defined C I = the number of additional angles subject to this correction, C and I>0 requires an additional line, format=(a31,20i2), specifying C the numbers of atoms carrying such angles C V = optional initial value for this parameter C C N = -8 isotopic bond length change relative to the value in C the parent - this is only for use with Deuterium C substitution. C M = the atom for which this bond is defined C I = the number of additional bonds subject to this correction, C and I>0 requires an additional line, format=(a31,20i2), specifying C the numbers of atoms carrying such bonds C V = optional initial value for this parameter C C FIX = 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. C This is one of only two cases for which the text in C in columns 1-31 is relevant (the other is XXX for an excluded C spectroscopic constant). Otherwise this field is ignored and C can be used for arbitrary explanatory comments. C C C------------------------------------------------------------------------------- C lines17-23 - block defining the spectroscopic constants to be fitted to C------------------------------------------------------------------------------- C C line 17: The total number of available spectroscopic constants. C C line 18+: DEFINITION OF THE SPECTROSCOPIC CONSTANT TO BE FITTED TO: C N,M,B,dI,dB C C read in with format=(31x,2i2,3f14.0) so that the numbers have C to be EXACTLY in the columns indicated below: C C column: 32 36 50 64 78 C | | | | | C ++--++++++++++++++--------------++++++++++++++ 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 Constant,species,value....XXX 1 1 312142.390 NH2CN C Constant,species,value....... 1 1 312142.390 0.1 30000. NH2CN C C N - the code defining the type of spectroscopic 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 where P.a=(1/2)[Ib+Ic-Ia] etc. C C 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 M = the number of the isotopic species C B = the value Bobs of the experimentally determined constant, in MHz for C rotational constants and uA**2 for planar moments. C The constants should be declared in increasing order of C isotopic species number, although all constants for C a given species are not required. c c It is optionally possible to specify EITHER fourth OR fifth numbers, c dI or dB (usually only one allowed), such that: C C dI - vibration-rotation correction dI=Icorr-Iobs (in uA**2) C to the experimental moment of inertia defined by N and B (even if C experimental rotational constants are specified). In this case C the dB field should be left blank. C C dB - vibration-rotation correction dB=Bcorr-Bobs (in MHz) C to the experimental rotational constant defined by N and M. C C --> this is most often (Be-B0) from anharmonic force field calculations C and used for r_e^SE structure evaluation C C In this case the space corresponding to dI is normally to be left blank. C BUT if a value is placed in this position then it is assumed to be C the (dimensionless) rotational g-factor, g_bb, for rotational constant B. C STRFIT will then subtract from the pertinent rotational constant C an electronic correction dB_el (in MHz) = 0.000544617 * g_bb * B, C where B is taken to be B+dB. C C The structure will then be fitted to corrected quantities C Icorr=Iobs+dI, Iobs=conv/Bobs or to C Icorr=conv/(Bobs+dB) or to C Icorr=conv/(Bobs+dB-dB_el) C where Iobs is the moment of inertia derived from the preceding C rotational constant, or to P=Pobs+dI if planar moment has been C specified C C The main purpose of the dI option is evaluation of r* (r_z) C structures by using dI's from a harmonic force field C C The main purpose of the dB option is evaluation of C semi-experimental equilibrium (r_e^SE) structures to either C the specified B = B0_exp + (Be-B0)_calc C or to B = B0_exp + (Be-B0)_calc - dB_el C C XXX = Any of the input rotational constants can be excluded from the fit C by writing XXX in columns 29-31 of the declaration line. C This is one of two only two cases for which the text in C in columns 1-31 is relevant. Otherwise this field is ignored and C can be used for arbitrary explanatory comments. C C C------------------------------------------------------------------------------- C lines24-26 - Block defining the isotopic changes C------------------------------------------------------------------------------- C C There is to be one such block for each substituted isotopic species C represented in the spectroscopic constants. Species 1 is taken to C be the parent species, the first isotopic block is to be for C species number 2, and successive such blocks are to be for C successively numbered species. C Each block begins with a repeat count line (line 24), which is C followed by the declaring lines. C C line 24: The number of changes relative to the parent for a given C isotopic species C C line 25+ - DEFINITION OF ISOTOPIC CHANGE: C N,M,V C C read in with format=(31x,2i2,3f14.0) so that the numbers have C to be EXACTLY in the columns indicated below: C C column: 32 36 49 C | | | C ++--++++++++++++++ C C N = the number of the isotopic atom C M = the number of the changed parameter (1,2,3,4 for C bond length, bond angle, dihedral angle and atomic mass resp.), C V = the value of this parameter in the substituted species. C When the parameter is not a mass (N not equal to 4) then V declares C the change relative to the value in the parent species C C------------------------------------------------------------------------------- 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 1-31 can be used for any clarifying comments, which are C up to the user and are not echoed to the output C C - Columns to the right of the input fields can be also be used for C comments, which are not echoed to the output. C Usually columns 78 and higher are safe for this purpose except in C the definitions of parameters of fit where the descriptive C alphanumeric parameter extends up to and inclusive of column 91 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 it C C - Additional blocks of comments are allowed before all lines C with a repeat count: i.e. 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 comment lines have to have a '!' character in the first C column, and there is no limit on their number. Such lines are C echoed at the top of the output file. 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 the resulting possibility of perpendicular C axis switching. Steps have been taken to stabilise such behaviour, C and those are marked by the string May2009 in the code C 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 the 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 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 data files from C the webpage 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 Command line example for gfortran: C gfortran -fno-automatic strfit.for -o strfit C C Command line example for Intel Visual Fortran 9: C ifort -nopdbfile -nodebug -traceback -arch=pn1 -tune=pn1 -O2 -Qsave C -ccdefault:fortran -fpscomp:filesfromcmd strfit.for C C Command line example for Intel Visual Fortran 11: C ifort -nopdbfile -nodebug -traceback -O3 -Qsave -ccdefault:fortran C -fpscomp:filesfromcmd strfit.for C C____________________________________________________________________________ C C IMPLICIT REAL(8) (A-H,O-Z) EXTERNAL LMUSER C parameter (maxpar=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,maxato=300,maxcor=maxato*3,naisot=35) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /INPUT/CARTRF(MAXATO,4), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2),NPARF(maxpar,3), * NREP(maxpar,maxrep),NISOT,NPAR,NOBSC COMMON /exclc/VCEXCL(maxcon,4),NCEXCL(maxcon,2),NEXCL COMMON /CHRDAT/NAMPAR,DESPAR,COMENT COMMON /DERIV/DER(MAXCON,maxpar) COMMON /coords/CORDP(MAXATO,3) COMMON /COV/COVTMP(maxpar,maxpar) COMMON /EPSIL/ristr(maxcon,3),ritot(maxcon,3) COMMON /files/filnam,filstf,filout,filpmi,filxyz c CHARACTER(100) filnam,filstf,filout,filpmi,filxyz 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,despar(maxpar)*20 CHARACTER(2) DTNAM(NAISOT),atnam REAL(8) DELTA(9),CORDP0(MAXATO,3),DERPC(maxcor,maxpar), * CCMAT(maxpar,maxpar), * ZT(maxcor,maxpar),ZTZP(maxcor,maxcor),topicm(naisot) 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 DATA nameco/' a ',' b ',' c ','b+c','a+b','P.a','P.b','P.c'/ DATA TOPICM * /1.0078250321D0, 2.014101778D0, 3.016049268D0, 12.0000000D0, ! H H H C * 13.0033548378D0, 28.976496D0, 14.0030744D0, 15.0001077D0, ! C Si N N * 15.9949150D0, 16.999133D0, 17.999160D0, 18.9984046D0, ! O O O F * 31.97207100D0, 33.96786690D0, 34.96885268D0, 36.96590259D0, ! S S Cl Cl * 78.918329D0, 80.916292D0, 126.904473D0, 6.01513D0, ! Br Br I Li * 7.01601D0, 9.01219D0, 10.01294D0, 11.009305D0, ! Li Be B B * 27.97693D0, 30.973765D0, 19.9924405D0, 21.9913847D0, ! Si P Ne Ne * 35.9675445D0, 39.9623842D0, ! Ar Ar * 81.913482D0, 82.9141314D0, 83.9115034D0, 85.9106159D0, ! Kr Kr Kr Kr * 1.0D-10/ ! * DATA DTNAM * /' H', ' H', ' H', ' C', * ' C', 'Si', ' N', ' N', * ' O', ' O', ' O', ' F', * ' S', ' S', 'Cl', 'Cl', * 'Br', 'Br', ' I', 'Li', * 'Li', 'Be', ' B', ' B', * 'Si', ' P', 'Ne', 'Ne', * 'Ar', 'Ar', * 'Kr', 'Kr', 'Kr', 'Kr', * ' *'/ 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',I3,': Chi-squared =',F17.8, 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 NITFIN=niter ALAFIN=ALAMDA 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 MOLNAM.PMI C (with masses of the parent isotopologue) C 3/ to UNIT=2 echo of fitted structure to MOLNAM.XYZ 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=filpmi(1:len_trim(filpmi)),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 1057 j=1,3 ! Fortran2018 cartdf(i,j)=cartrf(i,j) 1057 continue ! Fortran2018 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) c c___standard repetitions IF(NPARF(NN,3).GT.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 c c___repetitions with changes IF(NPARF(NN,3).LT.0)THEN DO 522 NNN=1,-NPARF(NN,3) if(nrep(nn,nnn).gt.0)then CARTDF(NREP(NN,NNN),NPARF(NN,2))= * VPARF(NN)+vdiff(nn,nnn) else CARTDF(-NREP(NN,NNN),NPARF(NN,2))= * -VPARF(NN)+vdiff(nn,nnn) endif 522 CONTINUE ENDIF c endif 120 CONTINUE C C...write current CART data to MOLNAM.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 lc=0 do 6104 k=1,npar if(len_trim(despar(k)).gt.lc)lc=len_trim(despar(k)) 6104 continue c if(ALAFIN.lt.1.d-7)then write(3,'(/1X,78(1H_)/1x,a,i4,a,E9.2/)')'fit after:',NITFIN, 1 ' iterations, ALAMDA= ',ALAFIN else write(3,'(/1X,78(1H_)/1x,a,i4,a,E9.2,a/)')'fit after:',NITFIN, 1 ' iterations, ALAMDA= ',ALAFIN, 1 ' <-- WARNING: poor convergence!' endif c WRITE(3,'('' 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 ! fitted parameter n=n+1 ERRPAR(N)=SQRT(COVAR(N,N))*SWDEV WRITE(*,6032)NAMPAR(k),A(N),ERRPAR(N),despar(k)(1:lc) c if(nparf(k,3).eq.0)then WRITE(3,6032)NAMPAR(k),A(N),ERRPAR(N),despar(k)(1:lc) WRITE(4,6032)NAMPAR(k),A(N),ERRPAR(N),despar(k)(1:lc) endif c c___standard repetitions if(nparf(k,3).gt.0)then WRITE(3,6033)NAMPAR(k),A(N),ERRPAR(N), * despar(k)(1:lc), * (nrep(k,nn),nn=1,nparf(k,3)) WRITE(4,6033)NAMPAR(k),A(N),ERRPAR(N), * despar(k)(1:lc), * (nrep(k,nn),nn=1,nparf(k,3)) endif c c___repetitions with changes if(nparf(k,3).lt.0)then WRITE(3,6033)NAMPAR(k),A(N),ERRPAR(N), * despar(k)(1:lc), * (nrep(k,nn),nn=1,-nparf(k,3)) WRITE(4,6033)NAMPAR(k),A(N),ERRPAR(N), * despar(k)(1:lc), * (nrep(k,nn),nn=1,-nparf(k,3)) endif c else ! fixed parameter if(nparf(k,3).eq.0)then WRITE(3,6132)NAMPAR(k),vparf(k),despar(k)(1:lc) WRITE(4,6132)NAMPAR(k),vparf(k),despar(k)(1:lc) endif c c___standard repetitions if(nparf(k,3).gt.0)then WRITE(3,6133)NAMPAR(k),vparf(k), * despar(k)(1:lc), * (nrep(k,nn),nn=1,nparf(k,3)) WRITE(4,6133)NAMPAR(k),vparf(k), * despar(k)(1:lc), * (nrep(k,nn),nn=1,nparf(k,3)) endif c c___repetitions with changes if(nparf(k,3).lt.0)then WRITE(3,6133)NAMPAR(k),vparf(k), * despar(k)(1:lc), * (nrep(k,nn),nn=1,-nparf(k,3)) WRITE(4,6133)NAMPAR(k),vparf(k), * despar(k)(1:lc), * (nrep(k,nn),nn=1,-nparf(k,3)) endif c endif 6004 CONTINUE 6032 FORMAT(1X,A,F13.6,' +-',F10.6,3x,a) 6033 FORMAT(1X,A,F13.6,' +-',F10.6,3x,a,' and at atom',20I3) 6132 FORMAT(1X,A,' [',F11.6,' ] FIXED ',3x,a) 6133 FORMAT(1X,A,' [',F11.6,' ] FIXED ',3x,a, * ' and at atom',20I3) C c...Overall statistics c chisqmhz=0.d0 do 40 i=1,NPT chisqmhz=chisqmhz+(vcobs(i,3)-vcobs(i,4))**2 40 continue swdevmhz=dsqrt( chisqmhz/ndegf ) c WRITE(3,6005)CHISQ,SWDEV,ndegf,swdevmhz WRITE(4,6005)CHISQ,SWDEV,ndegf,swdevmhz WRITE(*,6005)CHISQ,SWDEV,ndegf,swdevmhz WRITE(*,'(1X/)') 6005 FORMAT(/' Chi-squared =',F20.10, * ' = Sum( (Iobs-calc)**2 )'/ * ' Deviation of fit =',F16.6, * ' uA^2 = Sqrt(Chisq/Ndegf), Ndegf=',i3// * ' Note that the fit is to moments of inertia but it also ', * 'corresponds to:'/ * ' Deviation of fit =',F16.6, * ' MHz = Sqrt( Sum( (Bo-c)**2 )/Ndegf ) ') 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',8X,'Bcalc',7X,'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,4X,2F13.5,F11.5) C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C Obs. and Obs.-Calc. quantities for constants excluded from fit C C (the final fit has already set up the updated CARTDF and RMC,RMD C so that only isotopic changes in the masses are required) C if(nexcl.gt.0)then WRITE(3,5) WRITE(3,23) WRITE(4,5) WRITE(4,23) 5 FORMAT(1x/' Deviations for constants excluded from the fit:') DO 220 I=1,NEXCL c DO 3 NN=1,NATMS CARTDF(NN,4)=CARTRF(NN,4) 3 CONTINUE niso=ncexcl(i,2) do 4 nn=1,numchg(niso) if(nchang(niso,nn,2).eq.4)then cartdf(nchang(niso,nn,1),4)=vchang(niso,nn) endif 4 continue c call momin(ncexcl(i,2),0) ! <----- vcexcl(i,2)=PMI(ncexcl(i,1),1) vcexcl(i,4)=PMI(ncexcl(i,1),2) WRITE(3,22)ncexcl(i,2),nameco(ncexcl(i,1)), 1 VCEXCL(I,1),VCEXCL(I,2),VCEXCL(I,1)-VCEXCL(I,2), 1 VCEXCL(I,3),VCEXCL(I,4),VCEXCL(I,3)-VCEXCL(I,4) 220 CONTINUE endif C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 CCMAT(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),(CCMAT(N,L),L=M,LLLL) WRITE(4,303)n,cnames(N)(1:19),(CCMAT(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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C Final Cartesians and their estimated uncertainties written C to main output C C...zeroth order run, CORDP0 used to preserve the reference principals C do 300 i=1,natms cartdf(i,4)=cartrf(i,4) 300 continue call momin(1,0) ! <------- do 310 i=1,natms do 1310 j=1,3 ! Fortran2018 cordp0(i,j)=cordp(i,j) 1310 continue ! Fortran2018 310 continue c c...evaluate the derivatives to serve as the Z matrix: C C N - the number of the parameter of fit among fitted parameters C NV - the number of parameter of fit among declared parameters C I - the number of atom c II - the number of the Cartesian coordinate, in order a1,b1,c1,a2,b2,c2,a3... c do 313 i=1,3*natms do 1313 j=1,NPAR ! Fortran2018 derpc(i,j)=0.d0 1313 continue ! Fortran2018 313 continue do 330 i=1,ninfit do 1330 j=1,ninfit ! Fortran2018 covar(i,j)=swdsq*covar(i,j) 1330 continue ! Fortran2018 330 continue c n=0 do 311 nv=1,NPAR if(ifit(nv).eq.0)goto 311 n=n+1 c c...increment the required internal taking care of multiplicity and C possible symmetrically opposed declaration c if(nparf(nv,1).gt.0)then iatm=nparf(nv,1) iint=nparf(nv,2) cartdf(iatm,iint)=cartdf(iatm,iint)+delta(iint) c c___standard repetitions IF(NPARF(NV,3).GT.0)THEN DO 326 NNN=1,NPARF(NV,3) I=iabs(NREP(NV,NNN)) CARTDF(I,iint)=CARTDF(I,iint)+DELTA(iint) if(nrep(nv,nnn).lt.0) * cartdf(i,iint)=cartdf(i,iint)-2*delta(iint) 326 CONTINUE ENDIF c c___repetitions with changes IF(NPARF(NV,3).LT.0)THEN DO 526 NNN=1,-NPARF(NV,3) I=iabs(NREP(NV,NNN)) CARTDF(I,iint)=CARTDF(I,iint)+DELTA(iint) if(nrep(nv,nnn).lt.0) * cartdf(i,iint)=cartdf(i,iint)-2*delta(iint) 526 CONTINUE ENDIF c endif c c...increment a non-internal c if(nparf(nv,1).lt.0)then endif c c...principals for incremented internal c CALL momin(1,0) ! <------- c c...decrement the internal c if(nparf(nv,1).gt.0)then cartdf(iatm,iint)=cartdf(iatm,iint)-delta(iint) c c___standard repetitions IF(NPARF(NV,3).GT.0)THEN DO 327 NNN=1,NPARF(NV,3) I=iabs(NREP(NV,NNN)) CARTDF(I,iint)=CARTDF(I,iint)-DELTA(iint) if(nrep(nv,nnn).lt.0) * cartdf(i,iint)=cartdf(i,iint)+2*delta(iint) 327 CONTINUE ENDIF c c___repetitions with changes IF(NPARF(NV,3).LT.0)THEN DO 5327 NNN=1,-NPARF(NV,3) I=iabs(NREP(NV,NNN)) CARTDF(I,iint)=CARTDF(I,iint)-DELTA(iint) if(nrep(nv,nnn).lt.0) * cartdf(i,iint)=cartdf(i,iint)+2*delta(iint) 5327 CONTINUE ENDIF c endif c ii=0 do 312 i=1,natms DO 1312 j=1,3 ! Fortran2018 ii=ii+1 if(nparf(nv,1).gt.0)then derpc(ii,n)=(cordp(i,j)-cordp0(i,j))/delta(iint) endif 1312 continue ! Fortran2018 312 continue c 311 continue c c...diagnostic output c c write(3,'(1x/a/)') c * ' Last two columns are CHECK=SUMj(Der(i,j)*par(j)) '// c * 'and the relevant Cartesian' c do 314 ii=1,3*natms c iatm=(ii-1)/3+1 c icrd=ii-(iatm-1)*3 c check=0.d0 c n=0 c do 315 nv=1,npar c if(ifit(nv).eq.0)goto 315 c n=n+1 c check=check+derpc(ii,n)*a(n) c315 continue c write(3,'(50(F12.6))')(derpc(ii,n),n=1,ninfit), c * check,cordp0(iatm,icrd) c314 continue c c c...Evaluate the variance covariance matrix of the coordinates, eq.46 on c p.56 of Albritton,Schmeltekopf,Zare, "An introduction to the least- c squares fitting of spectroscopic data", in Molecular Spectroscopy: c Modern research, Vol.II, K.N.Rao (Ed.), Academic Press 1976. C C Z x THETA x Z' , which in this case is C C DERPC x COVAR*SWDSQ x DERPC' C C dimensions = (3*natms x 3*natms) C = (3*natms x ninfit) x (ninfit x ninfit) * (ninfit x 3*natms) C c Note that the array COVAR out of MRQMIN is the dispersion matrix V of c eq.11,12 p.17 ASZ, and is related to the variance-covariance matrix THETA c of the parameters by THETA=sigma^2*COVAR c do 316 i=1,3*natms do 1316 j=1,ninfit ! Fortran2018 sum=0.d0 do 317 ii=1,ninfit sum=sum+derpc(i,ii)*covar(ii,j) 317 continue zt(i,j)=sum 1316 continue ! Fortran2018 316 continue c do 321 i=1,3*natms do 1321 j=1,3*natms ! Fortran2018 sum=0.d0 do 322 ii=1,ninfit sum=sum+zt(i,ii)*derpc(j,ii) 322 continue ztzp(i,j)=sum 1321 continue ! Fortran2018 321 continue c C c...Final principal coordinates c 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,CORDp0(I,1),CORDp0(I,2),CORDp0(I,3),cartrf(I,4) 3112 CONTINUE 3114 FORMAT(2X,I2,2X,3F15.6,F17.7) C c...and their errors c write(3,320) 320 FORMAT(1X/' Principal coordinates and estimated uncertainties:'// * ' ATOM NO. A dA B dB', * ' C dC MASS'/) do 318 i=1,natms ii=(i-1)*3 write(3,319)i, * cordp0(i,1),dsqrt(ztzp(ii+1,ii+1)), * cordp0(i,2),dsqrt(ztzp(ii+2,ii+2)), * cordp0(i,3),dsqrt(ztzp(ii+3,ii+3)),cartrf(I,4) 318 continue 319 format(i4,1x,3(F13.5,f9.5),f16.6) write(3,325) 325 format(' | |', *' |'/ *' NOTES: 1/ only the uncertainties for those coordinates ', *'which are completely'/ *' defined by the fitted internals should be trusted'/ *' 2/ the uncertainties are somewhat limited by the ', *'linear approximation'/ *' coord=(d coord/d parameter)*parameter used for ', *'evaluation'/ *' 3/ only the effect of the internals R, A, and D ', *'is propagated'/) WRITE(3,'(1X,78(1H_))') C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...Final inertial contributions from c,d,Laurie and shift parameters C (run an additional calculation without non-internals and subtract C from the total) 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 isotopologue specifications to MOLNAM.PMI c IF(NISOT.gt.1)then DO 6007 N=2,NISOT WRITE(4,6042)N DO 7007 NN=1,NUMCHG(N) ! Fortran2018 WRITE(4,6043)' atom no.,parameter no.,value ', 1 (NCHANG(N,NN,I),I=1,2),VCHANG(N,NN) 7007 CONTINUE ! Fortran2018 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...write the MOLNAM.XYZ file C open(2,file=filxyz(1:len_trim(filxyz)),status='unknown',err=521) c WRITE(2,'(i5)')NATMS WRITE(2,'(A)')COMENT DO 3208 N=1,NATMS do 3212 nn=1,naisot if(nint(cartrf(n,4)).eq.nint(topicm(nn)))goto 3214 3212 continue 3214 atnam=dtnam(nn) if(atnam.eq.' *')then atnam='He' write(*,'(1x,a,i4,3a)')'NOTE: dummy atom',n, * ' converted to He in ',filxyz(1:len_trim(filxyz)),' output' endif WRITE(2,3211)atnam,CORDp0(n,1),CORDp0(n,2),CORDp0(n,3) 3208 CONTINUE 3211 format(a2,3F15.6) c close(2) c write(*,'(1x//'' Full output written to: '',a)') * filout(1:len_trim(filout)) write(*,'( '' Output for PMIFST written to: '',a)') * filpmi(1:len_trim(filpmi)) write(*,'( '' .XYZ output written to: '',a)') * filxyz(1:len_trim(filxyz)) c STOP c c 520 write(*,'(1X//'' **** ERROR on opening file: '',a//)') * filpmi(1:len_trim(filpmi)) STOP c 521 write(*,'(1X//'' **** ERROR on opening file: '',a//)') * filxyz(1:len_trim(filxyz)) STOP c 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 connection 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-parameter 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 10 for d_a = 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 VDIFF - differences in parameter values for additional occurrences of C a given parameter, with indexing as for NREP 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 NEXCL - the number of declared constants to be excluded from fit C NCEXCL - same as NCOBS but for excluded constants C VCEXCL - same as VCOBS but for excluded constants C 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 - automatic alphanumeric descriptor for the parameter of fit, C set up in DATAIN and used later C DESPAR - optional, user supplied additional alphanumeric descriptor for c the parameter of fit C............................................................................ C IMPLICIT REAL(8) (A-H,O-Z) C parameter (maxpar=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,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), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2),NPARF(maxpar,3), * NREP(maxpar,maxrep),NISOT,NPAR,NOBSC COMMON /exclc/VCEXCL(maxcon,4),NCEXCL(maxcon,2),NEXCL COMMON /CHRDAT/NAMPAR,DESPAR,COMENT COMMON /files/filnam,filstf,filout,filpmi,filxyz CHARACTER CNAME(8)*3,PARNM(6)*2,dastr*26 CHARACTER(20) NAMPAR(maxpar),despar(maxpar) CHARACTER COMENT*78,line*200,text*31,isonam*3,endcom*10 CHARACTER(100) filnam,filstf,filout,filpmi,filxyz,errmes 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 line and the CART definition scheme for the molecule - this C defines the parent species C WRITE(*,'(1x//)') WRITE(*,6002) 6002 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | STRFIT - General structure fitting', * ' program using CART definitions',T79,'|'/ * ' |',76(1H_),'|'/' version 14.VI.2021',T64,'Zbigniew KISIEL'//) C write(*,2)'CURRENT LIMITS:' write(*,2)'' write(*,2)'different isotopic species =',maxiso write(*,2)'isotopic changes relative to the parent =',maxchg write(*,2)'structural parameters of fit =',maxpar write(*,2)'repetitions of a given structural parameter =',maxrep write(*,2)'spectroscopic constants to fit to =',maxcon write(*,2)'atoms in the structure =',maxato write(*,2)'' 2 format(29x,a,i5) C 6001 WRITE(*, 3040,advance='NO') * 'Name of data file (implicit extension .STF): ' 3040 FORMAT(1X/1X,A) READ(*,'(A)',ERR=6001)FILNAM filstf=filnam(1:len_trim(filnam))//'.stf' OPEN(2,FILE=FILSTF(1:len_trim(filstf)),ERR=6001,status='old') nlin=0 C filout=filnam(1:len_trim(filnam))//'.out' filpmi=filnam(1:len_trim(filnam))//'.pmi' filxyz=filnam(1:len_trim(filnam))//'.xyz' c OPEN(3,FILE=FILOUT(1:len_trim(filout)),STATUS='UNKNOWN',ERR=6010) C errmes='reading the comment' READ(2,'(A)',ERR=6081,end=6080)LINE nlin=nlin+1 COMENT=LINE(1:len_trim(line)) write(3,6002) WRITE(3,6703)COMENT 6703 FORMAT(1x,78(1H_)//1X,A/1X,78(1H_)/) C errmes='skipping comment lines' 6702 READ(2,'(a)',err=6081,end=6080)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE(1:len_trim(line)) goto 6702 endif write(*,'(1x)') C c C------------------------------------------------------------------------------ C Read the CART deck C------------------------------------------------------------------------------ C errmes='reading the number of atoms' READ(LINE,3003,ERR=6081,end=6080)NATMS,idump 3003 FORMAT(2I5) DO 3206 N=1,NATMS errmes='reading a CART atom definition' READ(2,'(A)',ERR=6081,end=6080)LINE nlin=nlin+1 READ(LINE,*,ERR=6081,end=6081) * (NCART(N,NN),NN=1,4),(CARTRF(N,NN),NN=1,4) c if(cartrf(n,4).lt.0.d0)then errmes='CART definition - negative value of atomic mass?' goto 6081 endif c if(ncart(n,2).lt.0.or.ncart(n,3).lt.0.or.ncart(n,4).lt.0)then errmes='CART atom numbers since all should be positive' goto 6081 endif c if(n.eq.1.and.(ncart(n,2)+ncart(n,3)+ncart(n,4)).ne.0)then errmes='CART atom 1 since all defining atoms have to be 0' goto 6081 endif c if(n.eq.2.and.(ncart(n,3)+ncart(n,4)).ne.0)then errmes='CART atom 2 since atoms NB and NC have to be 0' goto 6081 endif c if(n.eq.3.and.ncart(n,4).ne.0)then errmes='CART atom 3 since atom NC has to be 0' goto 6081 endif c if(ncart(n,1).ne.n)then errmes= * 'first CART number (to be the same as no of the defined atom)' goto 6081 endif c if(ncart(n,2).ge.n.or.ncart(n,3).ge.n.or.ncart(n,4).ge.n)then errmes= * 'CART atoms numbers NA,NB or NC (to be smaller than N)' goto 6081 endif c if( n.ge.4 .and. * ((ncart(n,2).eq.ncart(n,3)) .or. * (ncart(n,2).eq.ncart(n,4)) .or. * (ncart(n,3).eq.ncart(n,4))) )then errmes= * 'CART atoms since all have to be different' goto 6081 endif c 3206 CONTINUE c ndatms=0 DO 6020 N=1,NATMS if(abs(cartrf(n,4)).lt.0.00001d0)ndatms=ndatms+1 DO 7020 NN=1,4 ! Fortran2018 CARTDF(N,NN)=CARTRF(N,NN) 7020 CONTINUE ! Fortran2018 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------------------------------------------------------------------------------ C errmes='skipping comments and reading the number of parameters' 6700 READ(2,'(a)',ERR=6081,end=6080)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE(1:len_trim(line)) goto 6700 endif C NFIX=0 READ(line,6013,ERR=6081,end=6080)TEXT,NPAR ! number of parameters 6013 FORMAT(A,I10) WRITE(3,'(1x/1x,45(1h-))') WRITE(3,6011)'TOTAL NUMBER OF STRUCTURAL PARAMETERS: ',NPAR 6011 FORMAT(1X,A,I5) WRITE(3,'(1x,45(1h-))') WRITE(3,'(1x/1x,a)')'Parameters to be fitted: ' c DO 6004 N=1,NPAR errmes='reading parameter of fit definitions' READ(2,'(A)',ERR=6081,end=6080)LINE nlin=nlin+1 READ(LINE,5003,ERR=6081,end=6080)TEXT,(NPARF(N,I),I=1,3),value ! parameter definition 5003 FORMAT(A,3I2,f14.0) c if(nparf(n,1).eq.0)then errmes='parameter identifier which should be non-zero' goto 6081 endif c if(nparf(n,1).gt.natms.or.nparf(n,1).lt.-9)then errmes='parameter identifier which is out of range' goto 6081 endif c if(nparf(n,2).eq.0)then errmes= * 'parameter definition since zero value of second number' goto 6081 endif c if(nparf(n,1).gt.0.and.(nparf(n,2).gt.3.or.nparf(n,2).lt.0))then errmes= * 'definition of standard internal as second no should be 1,2,3' goto 6081 endif c if(nparf(n,1).eq.-4.and.(nparf(n,2).lt.1.or.nparf(n,2).gt.10)) * then errmes= * 'definition of c parameter as second number is out of range' goto 6081 endif c if( (nparf(n,1).eq.-6.or.nparf(n,1).eq.-7 * .or.nparf(n,1).eq.-9) . and. * (nparf(n,2).lt.1 .or.nparf(n,2).gt.natms))then errmes= * 'Laurie parameter, second number out of range' goto 6081 endif c if( nparf(n,1).eq.-8 .and. * (nparf(n,2).lt.1 .or.nparf(n,2).gt.natms))then errmes= * 'isotopic change, second number out of range' goto 6081 endif c if(nparf(n,1).eq.-5.and.(nparf(n,2).lt.1 .or. * nparf(n,2).gt.10 .or. * nparf(n,2).eq.4 .or. * nparf(n,2).eq.5 .or. * nparf(n,2).eq.6 ))then errmes= * 'definition of d parameter as second number is out of range' goto 6081 endif c if(text(29:31).eq.'FIX')then ifit(n)=0 NFIX=NFIX+1 else ifit(n)=1 endif vparf(n)=value ! starting value if(len_trim(line).gt.71)then ! parameter alphanumeric despar(n)=line(72:) else despar(n)=' ' endif 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 C___if parameter value is undefined then get the starting value from C the CART deck if(value.eq.0.d0)VPARF(N)=CARTDF(NPARF(N,1),NPARF(N,2)) c 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...standard repetitions C IF(NPARF(N,3).GT.0)THEN errmes='reading repetitions for a standard parameter of fit ' READ(2,'(a)',err=6081,end=6080)LINE nlin=nlin+1 READ(LINE,6047,err=6081,end=6080)TEXT, * (NREP(N,NN),NN=1,NPARF(N,3)) 6047 FORMAT(A,20I2) c do 6083 nn=1,nparf(n,3) if(nrep(n,nn).eq.0.or.iabs(nrep(n,nn)).gt.natms)then errmes='atom number for a repeated parameter of fit' goto 6081 endif 6083 continue c 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 in the case of Laurie correction 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 C C...repetitions with changes = value incremented by a specified amount C IF(NPARF(N,3).LT.0)THEN errmes= * 'reading repetitions for a parameter of fit with changes' do 504 nn=1,iabs(nparf(n,3)) READ(2,'(a)',err=6081,end=6080)LINE nlin=nlin+1 READ(LINE(32:),*,err=6081,end=6080)NREP(N,NN),vdiff(n,nn) if(ifit(n).eq.1)then if(nrep(n,nn).gt.0)then WRITE(3,505)NREP(N,NN),vparf(n)+vdiff(n,nn), * vdiff(n,nn) else WRITE(3,505)NREP(N,NN),-vparf(n)+vdiff(n,nn), * vdiff(n,nn) endif endif 504 continue 505 FORMAT(27X,', also at atom ',I4,' = ',F10.3, * ' (difference = ',F10.3,' )') c ENDIF c 6004 CONTINUE ! NPAR loop C IF(NFIX.GT.0)THEN WRITE(3,'(1x/1x,a)') * 'Fixed parameters: ' do 6070 nn=1,npar if(ifit(nn).eq.0)then WRITE(3,6041)NAMPAR(NN),VPARF(NN) c c___standard repetitions IF(NPARF(NN,3).GT.0) * WRITE(3,6046)NPARF(NN,3),(NREP(NN,NNN),NNN=1,NPARF(NN,3)) c c___repetitions with changes IF(NPARF(NN,3).LT.0)then do 506 n=1,iabs(nparf(nn,3)) if(nrep(nn,n).gt.0)then WRITE(3,505)NREP(NN,N),vparf(nn)+vdiff(nn,n), * vdiff(nn,n) else WRITE(3,505)NREP(NN,N),-vparf(nn)+vdiff(nn,n), * vdiff(nn,n) endif 506 continue endif c endif 6070 continue ENDIF WRITE(3,'(1x,45(1h-)/)') C C------------------------------------------------------------------------------ C Read spectroscopic constants: nisot,nconst,Bobs,deltaI,deltaB C------------------------------------------------------------------------------ C errmes='skipping comments and reading the number'// * ' of spectroscopic constants' 6701 READ(2,'(a)',err=6081,end=6080)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE(1:len_trim(line)) goto 6701 endif C NISOT=1 lastiso=1 READ(LINE,6013,ERR=6081,end=6081)TEXT,NOBSC WRITE(3,'(1x)') WRITE(3,'(1x,86(1h-))') WRITE(3,'(1x,a,i5)')'TOTAL NUMBER OF SPECTROSCOPIC CONSTANTS: ', * NOBSC N=0 NEXCL=0 iboth=0 c WRITE(3,'(1x/a/a)')' Isotopic B_expt Ib_expt '// * ' dI dB dB_el B_corr Ib_corr', * ' species (or g_bb)' WRITE(3,'(1x,86(1h-))') c DO 6005 NN=1,NOBSC errmes='reading spectroscopic constants' READ(2,'(A)',ERR=6081,end=6080)LINE nlin=nlin+1 c c___spectroscopic constant declared but excluded from the fit c if(line(29:31).eq.'XXX')then errmes='reading excluded spectroscopic constants' NEXCL=NEXCL+1 READ(LINE,6003,ERR=6081,end=6080)TEXT,(NCEXCL(NEXCL,I),I=1,2), * VCEXCL(NEXCL,3),DELTAI,DELTAB c if(vcexcl(nexcl,3).le.0.0d0)then errmes='spectroscopic constant value, which should be > 0' goto 6081 endif c if(ncexcl(nexcl,1).lt.1.or.ncexcl(nexcl,1).gt.8)then errmes= * 'parameter describing the spectroscopic constant to be fixed' goto 6081 endif c if(DELTAI.NE.0.D0.AND.DELTAB.NE.0.D0)IBOTH=1 c VCEXCL(NEXCL,1)=CONV/VCEXCL(NEXCL,3) ! Ib from B, etc. if(ncexcl(nexcl,1).eq.4)vcexcl(nexcl,1)=4.d0*vcexcl(nexcl,1) ! B+C if(ncexcl(nexcl,1).eq.5)vcexcl(nexcl,1)=4.d0*vcexcl(nexcl,1) ! A+B if(ncexcl(nexcl,1).ge.6)then ! Pa,Pb,Pc vcexcl(nexcl,1)=vcexcl(nexcl,3) vcexcl(nexcl,3)=conv/vcexcl(nexcl,1) endif c bexpt= vcexcl(nexcl,3) riexpt=vcexcl(nexcl,1) dbel=0.d0 endcom=' excluded' c c...new/existing isotopic species c IF(NCexcl(Nexcl,2).gt.lastiso)then if(NCEXCL(Nexcl,2).gt.NISOT)NISOT=NCexcl(Nexcl,2) lastiso=ncexcl(nexcl,2) write(isonam,'(i3)')NISOT else isonam=' ' endif c c...all isotopic species c IF(DELTAI.NE.0.D0.AND.DELTAB.EQ.0.0D0)THEN VCexcl(Nexcl,1)=VCexcl(Nexcl,1)+DELTAI VCexcl(Nexcl,3)=CONV/VCexcl(Nexcl,1) WRITE(3,7038)isonam,CNAME(NCEXCL(Nexcl,1)),bexpt,riexpt, * deltai,deltab,dbel, * VCexcl(Nexcl,3),VCexcl(Nexcl,1), * endcom(1:len_trim(endcom)) endif c IF(DELTAB.NE.0.0d0)THEN VCexcl(Nexcl,3)=VCexcl(Nexcl,3)+DELTAB ! + (Be-B0) VCexcl(Nexcl,1)=CONV/VCexcl(Nexcl,3) if(deltai.ne.0.d0)then ! - DeltaB_electr dbel=0.000544617d0*DELTAI*VCexcl(Nexcl,3) VCexcl(Nexcl,3)=VCexcl(Nexcl,3)-dbel VCexcl(Nexcl,1)=CONV/VCexcl(Nexcl,3) endif WRITE(3,7038)isonam,CNAME(NCEXCL(Nexcl,1)),bexpt,riexpt, * deltai,deltab,dbel, * VCexcl(Nexcl,3),VCexcl(Nexcl,1), * endcom(1:len_trim(endcom)) ENDIF c IF(DELTAI.eq.0.d0.AND.DELTAB.eq.0.0d0)THEN WRITE(3,7038)isonam,CNAME(NCEXCL(Nexcl,1)),bexpt,riexpt, * deltai,deltab,dbel, * VCexcl(Nexcl,3),VCexcl(Nexcl,1), * endcom(1:len_trim(endcom)) ENDIF c goto 6005 endif ! .eq.'XXX' c c___spectroscopic constant for use in fitting c n=n+1 READ(LINE,6003,ERR=6081,end=6080)TEXT,(NCOBS(N,I),I=1,2), * VCOBS(N,3),DELTAI,DELTAB 6003 FORMAT(A,2I2,3F14.4) c if(vcobs(n,3).le.0.0d0)then errmes='spectroscopic constant value, which should be > 0' goto 6081 endif c if(ncobs(n,1).lt.1.or.ncobs(n,1).gt.8)then errmes= * 'parameter describing the spectroscopic constant to be fitted' goto 6081 endif c if(DELTAB.ne.0.d0.and.NCOBS(N,1).gt.4)then write(*,'(1x//'' **** ERROR: The use of dB is only allowed '', * ''for A,B,C,B+C.''/13x,''You need to correct this line:''// * 1x,a//)')line(1:len_trim(line)) stop endif c if(DELTAI.NE.0.D0.AND.DELTAB.NE.0.D0)IBOTH=1 c c Ib+Ic is evaluated from (B+C) by assuming B=C=(B+C)/2 c Ia+Ib is evaluated from (A+B) by assuming A=B=(A+B)/2 c c Note that these relations are approximations with quality dependent on c the asymmetry parameter c c For Pa,Pb,Pc the fictitious rotational constants 50579.01/Pa etc. are evaluated c as placeholders c VCOBS(N,1)=CONV/VCOBS(N,3) ! Ib from B, etc. if(ncobs(n,1).eq.4)vcobs(n,1)=4.d0*vcobs(n,1) ! B+C if(ncobs(n,1).eq.5)vcobs(n,1)=4.d0*vcobs(n,1) ! A+B if(ncobs(n,1).ge.6)then ! Pa,Pb,Pc vcobs(n,1)=vcobs(n,3) vcobs(n,3)=conv/vcobs(n,1) endif bexpt=vcobs(n,3) riexpt=vcobs(n,1) dbel=0.d0 endcom='' C C...new/existing isotopic species C IF(NCOBS(N,2).gt.lastiso)then if(NCOBS(N,2).gt.NISOT)NISOT=NCOBS(N,2) lastiso=ncobs(n,2) write(isonam,'(i3)')NISOT else isonam=' ' endif c c...all isotopic species c IF(DELTAI.NE.0.D0.and.DELTAB.eq.0.d0)THEN VCOBS(N,1)=VCOBS(N,1)+DELTAI VCOBS(N,3)=CONV/VCOBS(N,1) WRITE(3,7038)isonam,CNAME(NCOBS(N,1)),bexpt,riexpt, * deltai,deltab,dbel, * VCOBS(N,3),VCOBS(N,1),endcom(1:len_trim(endcom)) 7038 FORMAT(a,2x,a,f13.5,f11.5,f10.5,f11.5,f9.5,f13.5,f11.5,a) ENDIF c IF(DELTAB.NE.0.0d0)THEN VCOBS(N,3)=VCOBS(N,3)+DELTAB ! + (Be-B0) VCOBS(N,1)=CONV/VCOBS(N,3) if(DELTAI.NE.0.d0)then dbel=0.000544617d0*DELTAI*VCOBS(N,3) VCOBS(N,3)=VCOBS(N,3)-dbel ! - DeltaB_electr VCOBS(N,1)=CONV/VCOBS(N,3) ENDIF if(ncobs(n,1).eq.4)vcobs(n,1)=4.d0*vcobs(n,1) WRITE(3,7038)isonam,CNAME(NCOBS(N,1)),bexpt,riexpt, * deltai,deltab,dbel, * VCOBS(N,3),VCOBS(N,1),endcom(1:len_trim(endcom)) ENDIF c IF(DELTAI.eq.0.d0.AND.DELTAB.eq.0.0d0)THEN WRITE(3,7038)isonam,CNAME(NCOBS(N,1)),bexpt,riexpt, * deltai,deltab,dbel, * VCOBS(N,3),VCOBS(N,1),endcom(1:len_trim(endcom)) ENDIF c 6005 CONTINUE c if(iboth.eq.1)then write(*,'(1x,a/1x,a)') * '***** NOTE ------> there are cases where both dI and dB '// * 'are non-zero:', * ' assuming dI is g_bb' endif c if(n.lt.nobsc)then WRITE(3,'(1x)') WRITE(3,6011)' NO OF CONSTANTS TO FIT TO: ',N WRITE(3,'(1X,A,I5)') * ' NO OF EXCLUDED CONSTANTS: ',nexcl NOBSC=N endif c WRITE(3,'(1x,86(1h-))') write(3,'(a/a/a/a/)') * ' B_corr = B_expt + dB - dB_el, Ib_corr=505379.01/B_corr', * ' g_bb is dimensionles, dB_el= = 0.000544617 * g_bb * B, '// * 'where B=B_expt+dB', * ' or', * ' Ib_corr = Ib_expt + dI, B_corr =505379.01/Ib_corr' C C------------------------------------------------------------------------------ C Read definitions of isotopically substituted species C------------------------------------------------------------------------------ C NUMCHG(1)=0 IF(NISOT.EQ.1)GOTO 6006 write(3,'(1x/1x,55(1h-)/a/1x,55(1h-))') * ' DEFINITIONS OF SUBSTITUTED ISOTOPIC SPECIES' c DO 6007 N=2,NISOT errmes='skipping comments and reading no of isotopic changes' 6704 READ(2,'(a)',err=6081,end=6080)LINE nlin=nlin+1 if(line(1:1).eq.'!')then write(3,'(1x,a)')LINE(1:len_trim(line)) goto 6704 endif READ(LINE,6013,ERR=6081,end=6081)TEXT,NUMCHG(N) WRITE(3,6042)N 6042 FORMAT(' ISOTOPIC SPECIES',I3,', changes from ', 1 'parent species:') c errmes='reading the isotopic changes' DO 6017 NN=1,NUMCHG(N) READ(2,'(A)',ERR=6081,end=6080)LINE nlin=nlin+1 READ(LINE,6093,ERR=6081,END=6081) * TEXT,(NCHANG(N,NN,I),I=1,2),VCHANG(N,NN) c if(nchang(n,nn,1).le.0.or. * nchang(n,nn,1).gt.natms.or. * nchang(n,nn,2).le.0.or. * nchang(n,nn,2).gt.4.or. * vchang(n,nn).eq.0.d0)then errmes='isotopic change definition' goto 6081 endif c 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) WRITE(3,'(1x,55(1h-)/)') GOTO 6006 C 6080 WRITE(*,'(1X//'' **** ERROR: unexpected end of file while''/ * '' '',a//)')errmes(1:len_trim(errmes)) STOP C 6081 WRITE(*,'(1X//'' **** ERROR in:''/ * '' '',a// * '' --- Please check line'',i5, * '' of the input file:''/'' |''/ * '' --> '',A,A//)')errmes(1:len_trim(errmes)), * nlin,LINE(1:len_trim(line)),CHAR(7) STOP C 6006 CLOSE(2) c RETURN C C 6010 write(*,'(1X//'' **** ERROR on opening the input file: '',a//)') * filstf(1:len_trim(filstf)) STOP 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=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,maxato=300) C COMMON /DERIV/DER(MAXCON,maxpar) COMMON /INPUT/CARTRF(MAXATO,4), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2),NPARF(maxpar,3), * NREP(maxpar,maxrep),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,DESPAR,COMENT REAL(8) X,A(NA),Y,DYDA(NA) CHARACTER(20) NAMPAR(maxpar),DESPAR(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=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,maxato=300) C COMMON /MOMI/NCART(MAXATO,4),CARTDF(MAXATO,4),PMI(8,2),NATMS, * NDATMS,idump,niter COMMON /INPUT/CARTRF(MAXATO,4), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2),NPARF(maxpar,3), * NREP(maxpar,maxrep),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,DESPAR,COMENT COMMON /DERIV/DER(MAXCON,maxpar) REAL(8) A(MA) CHARACTER(20) NAMPAR(maxpar),DESPAR(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 RMC - current rm(1) c parameters set up from either VPARF or A C RMD - current rm(2) d parameters set up from either VPARF or 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=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,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), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2),NPARF(maxpar,3), * NREP(maxpar,maxrep),NISOT,NPAR,NOBSC COMMON /CHRDAT/NAMPAR,DESPAR,COMENT COMMON /DERIV/DER(MAXCON,maxpar) COMMON /EPSIL/ristr(maxcon,3),ritot(maxcon,3) REAL(8) A(MA),rlaur(maxpar),rlaurr(maxpar,maxrep) CHARACTER(20) NAMPAR(maxpar),DESPAR(maxpar) CHARACTER(78) COMENT REAL(8) DELTA(9) 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 3000 NNN=1,4 ! Fortran2018 CARTDF(NN,NNN)=CARTRF(NN,NNN) 3000 CONTINUE ! Fortran2018 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 c if(ifit(nn).eq.0)then ! fixed internal CARTDF(NPARF(NN,1),NPARF(NN,2))=VPARF(NN) c c___standard repetitions IF(NPARF(NN,3).GT.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 c c___repetitions with changes IF(NPARF(NN,3).LT.0)THEN DO 522 NNN=1,-NPARF(NN,3) if(nrep(nn,nnn).gt.0)then CARTDF(NREP(NN,NNN),NPARF(NN,2))= * VPARF(NN)+VDIFF(NN,NNN) else CARTDF(-NREP(NN,NNN),NPARF(NN,2))= * -VPARF(NN)+VDIFF(NN,NNN) endif 522 CONTINUE ENDIF c else ! fitted internal kk=kk+1 CARTDF(NPARF(NN,1),NPARF(NN,2))=A(kk) c c___standard repetitions IF(NPARF(NN,3).GT.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 c c___repetitions with changes IF(NPARF(NN,3).LT.0)THEN DO 532 NNN=1,-NPARF(NN,3) if(nrep(nn,nnn).gt.0)then CARTDF(NREP(NN,NNN),NPARF(NN,2)) * =A(kk)+VDIFF(NN,NNN) else CARTDF(-NREP(NN,NNN),NPARF(NN,2)) * =-A(kk)+VDIFF(NN,NNN) endif 532 CONTINUE ENDIF c endif c 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 call setrm(A,MA,ncall) ! <----- C C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...implement changes relative to parent for the current isotopic species C IF(N.gt.1)then 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 ENDIF C C...evaluate current molecular mass C 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 bond 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 (Laurie bond correction) c c___standard repetitions if(nparf(nn,3).gt.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 c c___repetitions with changes if(nparf(nn,3).lt.0)then endif c endif C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c...Laurie correction defined in terms of projection 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 (Laurie projection angle correction) c c___standard repetitions if(nparf(nn,3).gt.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 c c___repetitions with changes if(nparf(nn,3).lt.0)then endif c endif 200 continue C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...Laurie type correction for a bond angle C WARNING: do not use - only for experimenting with new terms and C 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 (Laurie angle correction) c c___standard repetitions if(nparf(nn,3).gt.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 c c___repetitions with changes if(nparf(nn,3).lt.0)then endif c 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 (deuteration stretch) c c___standard repetitions IF(NPARF(NN,3).GT.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 c c___repetitions with changes IF(NPARF(NN,3).LT.0)THEN ENDIF c endif 235 continue C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C 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 C___standard repetitions if(nparf(nv,3).gt.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 c c___repetitions with changes if(nparf(nv,3).lt.0)then endif c 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 C___standard repetitions if(nparf(nv,3).gt.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 C C___repetitions with changes if(nparf(nv,3).lt.0)then endif c 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 C___standard repetitions if(nparf(nv,3).gt.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 C C___repetitions with changes if(nparf(nv,3).lt.0)then endif c 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 C___standard repetitions IF(NPARF(nv,3).GT.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 C C___repetitions with changes IF(NPARF(nv,3).LT.0)THEN ENDIF c goto 2 endif C C...gradient for a standard internal: bond, bond angle, dihedral angle C CARTDF(I,J)=CARTDF(I,J)+DELTA(J) c c___standard repetitions IF(NPARF(NV,3).GT.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___repetitions with changes IF(NPARF(NV,3).LT.0)THEN DO 523 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) 523 CONTINUE ENDIF c C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...calculate moments of inertia for this (=Nth) isotopologue 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 c c___standard repetitions if(nparf(nn,3).gt.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 c c___repetitions with changes if(nparf(nn,3).lt.0)then endif c 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 NC=NC+1 IF(NC.GT.NOBSC)GOTO 17 GOTO 5 C 1 CONTINUE C C 17 RETURN END C C_____________________________________________________________________________ C SUBROUTINE SETRM(A,MA,ncall) c C To set up the RMC and RMD arrays of c and d parameters C c NCALL = -1: just zero the c and d arrays c = 0 or higher: set up the arrays C IMPLICIT REAL(8) (A-H,O-Z) parameter (maxpar=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,maxato=300,maxcor=maxato*3) c COMMON /ANCIL/RMC(6),RMD(3),rmci(6),rmdi(3) COMMON /INPUT/CARTRF(MAXATO,4), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2),NPARF(maxpar,3), * NREP(maxpar,maxrep),NISOT,NPAR,NOBSC REAL(8) A(MA) 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)return c kk=0 do 100 nn=1,npar if(ifit(nn).eq.1)kk=kk+1 if(nparf(nn,1).eq.-4)then ! c of rm(1) i=nparf(nn,2) c if(i.le.6)then ! c_i or c_ij if(ifit(nn).eq.1)then rmc(i)=a(kk) else rmc(i)=vparf(nn) endif endif c if(i.eq.7)then ! c_a=c_b 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 ! c_a=c_c 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 ! c_b=c_c 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 ! c_a=c_b=c_c 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 ! d of rm(2) i=nparf(nn,2) c if(i.le.6)then ! d_i if(ifit(nn).eq.1)then rmd(i)=a(kk) else rmd(i)=vparf(nn) endif endif c if(i.eq.7)then ! d_a=d_b 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 ! d_a=d_c 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 ! d_b=d_c 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 ! d_a=d_b=d_c 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 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: C C PMI(i,1) - moments of inertia C PMI(i,2) - rotational constants C RMCI(i) - inertial contribution from c terms C RMDI(i) - inertial contribution from d terms C C where i=1,2,3,4,5,6,7,8 stands for a,b,c,b+c,a+b,Pa,Pb,Pc: 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=99,maxiso=99,maxcon=3*maxiso,maxrep=20, * maxchg=20,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), * VPARF(maxpar),VCHANG(maxiso,maxchg),VDIFF(maxpar,maxrep), * VCOBS(MAXCON,4),IFIT(maxcon),NUMCHG(maxiso),NPARF(maxpar,3), * NCHANG(maxiso,maxchg,2),NCOBS(MAXCON,2), * NREP(maxpar,maxrep),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 c IF(CARTDF(3,2))5555,118,5555 IF(CARTDF(3,2).ne.0.d0)goto 5555 ! Fortran2018 CS=-0.333333333333D0 ! old label 118 SS= 0.942809041582D0 GO TO 121 5555 CS=DCOS(CON*CARTDF(3,2)) SS=DSIN(CON*CARTDF(3,2)) NA=NCART(3,2) c121 IF(NA-N1)128,122,128 121 IF(NA-N1.ne.0)goto 128 ! Fortran2018 CORD(3,1)=CORD(NA,1)+CARTDF(3,1)*CS ! old label 122 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) c IF(CARTDF(I,2))133,132,133 IF(CARTDF(I,2).ne.0.d0)goto 133 ! Fortran2018 CS=-0.33333333333D0 ! old label 132 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 6556 M=1,3 CORD(NO,M)=CORD(NA,M) DO 7556 K=1,3 ! Fortran2018 CORD(NO,M)=CORD(NO,M)+TRANS(M,K)*PRC(K) 7556 CONTINUE ! Fortran2018 6556 CONTINUE ! Fortran2018 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 isotopologue 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 isotopologue 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 isotopologue 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 1415 NN=1,3 ! Fortran2018 EIGINV(N,NN)=EIGINV(N,NN)/WW 1415 CONTINUE ! Fortran2018 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 3994 nn=1,3 ! Fortran2018 AA(N,NN)=0.0D0 3994 CONTINUE ! Fortran2018 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 1996 NN=1,3 ! Fortran2018 ABACK(N,NN)=AA(N,NN) 1996 CONTINUE ! Fortran2018 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 2000 nn=1,3 ! Fortran2018 if(n.ne.nn)aa(n,nn)=0.d0 aar(n,nn)=0.d0 2000 continue ! Fortran2018 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 2001 nn=1,3 ! Fortran2018 do 3001 nnn=1,3 ! Fortran2018 aar(n,nn)=aar(n,nn)+eiginv(n,nnn)*eiginv(nn,nnn) * *aa(nnn,nnn) 3001 continue ! Fortran2018 2001 continue ! Fortran2018 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 1999 NN=1,3 ! Fortran2018 AA(N,NN)=ABACK(N,NN) 1999 CONTINUE ! Fortran2018 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 isotopologue then: C C 1/ Use initial centre of mass 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 (for the parent 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 4060 NN=1,3 ! Fortran2018 EIGINV(N,NN)=EIGINV(N,NN)/WW 4060 continue ! Fortran2018 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 (for the parent only the diagonal terms are applicable) 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 RANGE=1.0D-12 !old label 5 c IF(MV-1) 10,25,10 IF(MV-1.eq.0)goto 25 ! Fortran2018 IQ=-N ! old label 10 DO 1020 J=1,N ! Fortran2018 IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0.0d0 c IF(I-J) 20,15,20 IF(I-J.ne.0)goto 20 ! Fortran2018 R(IJ)=1.0d0 ! old label 15 20 CONTINUE 1020 CONTINUE ! Fortran2018 C C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) C 25 ANORM=0.0d0 DO 1035 I=1,N ! Fortran2018 DO 35 J=I,N c IF(I-J) 30,35,30 IF(I-J.eq.0)goto 35 ! Fortran2018 IA=I+(J*J-J)/2 ! old label 30 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE 1035 CONTINUE ! Fortran2018 c IF(ANORM) 165,165,40 IF(ANORM.le.0.d0)goto 165 ! Fortran2018 ANORM=1.41421356237d0*DSQRT(ANORM) ! old label 40 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 c 62 IF(DABS(A(LM))-THR) 130,65,65 IF(DABS(A(LM))-THR.lt.0.d0)goto 130 ! Fortran2018 old label 62 IND=1 ! old label 65 LL=L+LQ MM=M+MQ X=0.5d0*(A(LL)-A(MM)) Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) ! old label 68 c IF(X) 70,75,75 IF(X.ge.0.d0)goto 75 ! Fortran2018 Y=-Y ! old label 70 75 SINX=Y/DSQRT(2.0d0*(1.0d0+(DSQRT(1.0d0-Y*Y)))) SINX2=SINX*SINX COSX=DSQRT(1.0d0-SINX2) ! old label 78 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 c IF(I-L) 80,115,80 IF(I-L.eq.0)goto 115 ! Fortran2018 c 80 IF(I-M) 85,115,90 IF(I-M.ge.0)then ! Fortran2018 old label 80 if(I-M.gt.0)goto 90 ! Fortran2018 goto 115 ! Fortran2018 endif ! Fortran2018 IM=I+MQ ! old label 85 GO TO 95 90 IM=M+IQ c 95 IF(I-L) 100,105,105 95 IF(I-L.ge.0)goto 105 ! Fortran2018 IL=I+LQ ! old label 100 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 c 115 IF(MV-1) 120,125,120 115 IF(MV-1.eq.0)goto 125 ! Fortran2018 ILR=ILQ+I ! old label 120 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 c 130 IF(M-N) 135,140,135 130 IF(M-N.eq.0)goto 140 ! Fortran2018 M=M+1 ! old label 135 GO TO 60 C C TEST FOR L = SECOND FROM LAST COLUMN C c 140 IF(L-(N-1)) 145,150,145 140 IF(L-(N-1).eq.0)goto 150 ! Fortran2018 L=L+1 ! old label 145 GO TO 55 c 150 IF(IND-1) 160,155,160 150 IF(IND-1.ne.0)goto 160 ! Fortran2018 IND=0 ! old label 155 GO TO 50 C C COMPARE THRESHOLD WITH FINAL NORM C c 160 IF(THR-ANRMX) 165,165,45 160 IF(THR-ANRMX.gt.0.d0)goto 45 ! Fortran2018 C C SORT EIGENVALUES AND EIGENVECTORS C 165 IQ=-N DO 1185 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 c IF(A(LL)-A(MM)) 170,185,185 IF(A(LL)-A(MM).ge.0.d0)goto 185 ! Fortran2018 X=A(LL) ! old label 170 A(LL)=A(MM) A(MM)=X c IF(MV-1) 175,185,175 IF(MV-1.eq.0)goto 185 ! Fortran2018 DO 180 K=1,N ! old label 175 ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) R(IMR)=X 180 CONTINUE ! Fortran2018 185 CONTINUE 1185 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 1111 J=1,N ! Fortran2018 IF(ATOP.GE.DABS(A(I,J)))GOTO 111 ATOP=DABS(A(I,J)) 1111 CONTINUE ! Fortran2018 111 CONTINUE c IF(ATOP)109,109,113 IF(ATOP.gt.0.d0)goto 113 RETURN ! old label 109 C C CALCULATE THE STOPPING CRITERION -- DSTOP C 113 AVGF=(N*(N-1))*0.55D0 D=0.0D0 DO 114 JJ=2,N DO 1114 II=2,JJ ! Fortran2018 S=A(II-1,JJ)/ATOP D=S*S+D 1114 CONTINUE ! Fortran2018 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 1130 JCOL=2,N ! Fortran2018 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 A(JCOL,J)=S*T+C*U 124 continue ! Fortran2018 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 1128 II=2,JJ ! Fortran2018 S=A(II-1,JJ)/ATOP D=S*S+D 1128 continue ! Fortran2018 128 CONTINUE DSTOP=(1.0D-06)*D 129 THRSH=DSQRT(D/AVGF)*ATOP 130 CONTINUE 1130 CONTINUE ! Fortran2018 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 The routine set for nonlinear least-squares fitting with the Levenberg- C Marquardt method is adapted from Press,Flannery,Teukolsky,Vetterling, C Numerical Recipes, The Art of Scientific Computing, Cambridge University C Press (1986) C_____________________________________________________________________________ C SUBROUTINE MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT, 1 COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA) 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 containing 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-squared 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 progress 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=99, 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 WRITE(*,'(1x/1x,a//)')'MRQMIN: Improper LISTA permutation' ! ...STRMOD STOP ! ...STRMOD ENDIF 12 CONTINUE IF (KK.NE.(MA+1)) THEN WRITE(*,'(1x/1x,a//)')'MRQMIN: Improper LISTA permutation' ! ...STRMOD STOP ! ...STRMOD ENDIF 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 1020 IZB=1,MA ! ...STRMOD Fortran2018 COVTMP(IZA,IZB)=COVAR(IZA,IZB) ! ...STRMOD 1020 continue ! Fortran2018 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=99, 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 significant c modification 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 WRITE(*,'(1x/1x,a//)')'GAUSSJ: Singular matrix' ! ...STRMOD STOP ! ...STRMOD 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) THEN WRITE(*,'(1x/1x,a//)') 'GAUSSJ: Singular matrix.' ! ...STRMOD STOP ! ...STRMOD ENDIF 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_____________________________________________________________________________