From f8314098f60219a0bcd9ad5a45b9910e0342b706 Mon Sep 17 00:00:00 2001 From: David Williams Date: Mon, 21 Oct 2024 14:01:28 +0200 Subject: [PATCH] initial commit. --- .gitignore | 3 + GNUmakefile | 13 +++ JTmod.incl | 38 ++++++++ cart2int.f | 263 ++++++++++++++++++++++++++++++++++++++++++++++++++ ctrans.f | 88 +++++++++++++++++ nnparams.incl | 188 ++++++++++++++++++++++++++++++++++++ tab2genetic.f | 48 +++++++++ 7 files changed, 641 insertions(+) create mode 100644 .gitignore create mode 100644 GNUmakefile create mode 100644 JTmod.incl create mode 100644 cart2int.f create mode 100644 ctrans.f create mode 100644 nnparams.incl create mode 100644 tab2genetic.f diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..1cf0d0a --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.o +/cart2int +/tab2genetic diff --git a/GNUmakefile b/GNUmakefile new file mode 100644 index 0000000..d8bb04c --- /dev/null +++ b/GNUmakefile @@ -0,0 +1,13 @@ +.SUFFIXES: +.SUFFIXES: .o .f + +FC=gfortran +FFLAGS=--check=bounds -Wall -fmax-errors=5 -Werror + +tab2genetic: tab2genetic.o ctrans.o cart2int.o + $(FC) $(FFLAGS) $^ -o $@ + +%.o : JTmod.incl nnparams.incl + +clean: + -rm -fr *.o diff --git a/JTmod.incl b/JTmod.incl new file mode 100644 index 0000000..277fbd5 --- /dev/null +++ b/JTmod.incl @@ -0,0 +1,38 @@ +*** Relevant parameters for the analytic model +*** offsets: +*** offsets(1): morse equilibrium (N-H) +*** offsets(2): reference angle (H-N-H) +*** offsets(3): -- +*** pat_index: vector giving the position of the +*** various coordinates (see below) +*** ppars: polynomial parameters for tmcs +*** vcfs: coefficients for V expressions. +*** wzcfs: coefficients for W & Z expressions. +*** ifc: inverse factorials. + + integer matdim + parameter (matdim=5) ! matrix is (matdim)x(matdim) + + real*8 offsets(2) + integer pat_index(maxnin) + +! NH3 params + parameter (offsets=[1.0228710942d0,120.d0]) + +!########################################################################## +! coordinate order; the first #I number of coords are given to the +! ANN, where #I is the number of input neurons. The position i in +! pat_index corresponds to a coordinate, the value of pat_index(i) +! signifies its position. +! +! The vector is ordered as follows: +! a,xs,ys,xb,yb,b,rs**2,rb**2,b**2, +! es*eb, es**3, eb**3,es**2*eb, es*eb**2 +! ri**2 := xi**2+yi**2 = ei**2; ei := (xi,yi), i = s,b +! +! parts not supposed to be read by ANN are marked by ';' for your +! convenience. +!########################################################################## +! a,rs**2,rb**2,es*eb,es**3,eb**3,es**2*eb,es*eb**2,b**2 #I=9 (6D) + parameter (pat_index=[1,10,11,12,13,14,2,3,9,6,4,5,7,8]) +!************************************************************************** diff --git a/cart2int.f b/cart2int.f new file mode 100644 index 0000000..189e523 --- /dev/null +++ b/cart2int.f @@ -0,0 +1,263 @@ + subroutine cart2int(cart,qint) + implicit none +! This version merges both coordinate transformation routines into +! one. JTmod's sscales(2:3) are ignored. +! This is the first version to be compatible with one of my proper 6D fits +! Time-stamp: <2024-10-22 13:52:59 dwilliams> + +! Input (cartesian) +! cart(:,1): N +! cart(:,1+i): Hi +! Output +! qint(1): a mode +! qint(2:3): e stretching modes +! qint(4:5): e bending modes +! qint(6): umbrella +! Internal Variables +! no(1:3): NO distances 1-3 +! pat_in: temporary coordinates +! axis: main axis of NO3 + include 'nnparams.incl' + include 'JTmod.incl' + + + real*8 cart(3,4),qint(maxnin) + + real*8 no(3), r1, r2, r3 + real*8 v1(3), v2(3), v3(3) + real*8 n1(3), n2(3), n3(3), tr(3) + real*8 ortho(3) + real*8 pat_in(maxnin) + logical ignore_umbrella,dbg_umbrella + logical dbg_distances + +!.. Debugging parameters +!.. set umbrella to 0 + parameter (ignore_umbrella=.false.) +! parameter (ignore_umbrella=.true.) + +!.. break if umbrella is not 0 + parameter (dbg_umbrella=.false.) +! parameter (dbg_umbrella=.true.) + +!.. break for tiny distances + parameter (dbg_distances=.false.) +! parameter (dbg_distances=.true.) + + integer k + +!.. get N-O vectors and distances: + do k=1,3 + v1(k)=cart(k,2)-cart(k,1) + v2(k)=cart(k,3)-cart(k,1) + v3(k)=cart(k,4)-cart(k,1) + enddo + no(1)=norm(v1,3) + no(2)=norm(v2,3) + no(3)=norm(v3,3) + +!.. temporarily store displacements + do k=1,3 + pat_in(k)=no(k)-offsets(1) + enddo + + do k=1,3 + v1(k)=v1(k)/no(1) + v2(k)=v2(k)/no(2) + v3(k)=v3(k)/no(3) + enddo + +!.. compute three normal vectors for the ONO planes: + call xprod(n1,v1,v2) + call xprod(n2,v2,v3) + call xprod(n3,v3,v1) + + do k=1,3 + tr(k)=(n1(k)+n2(k)+n3(k))/3.d0 + enddo + r1=norm(tr,3) + do k=1,3 + tr(k)=tr(k)/r1 + enddo + +! rotate trisector + call rot_trisec(tr,v1,v2,v3) + +!.. determine trisector angle: + if (ignore_umbrella) then + pat_in(7)=0.0d0 + else + pat_in(7)=pi/2.0d0 - acos(scalar(v1,tr,3)) + pat_in(7)=sign(pat_in(7),cart(1,2)) + endif + +!.. molecule now lies in yz plane, compute projected ONO angles: + v1(1)=0.d0 + v2(1)=0.d0 + v3(1)=0.d0 + r1=norm(v1,3) + r2=norm(v2,3) + r3=norm(v3,3) + do k=2,3 + v1(k)=v1(k)/r1 + v2(k)=v2(k)/r2 + v3(k)=v3(k)/r3 + enddo + +! make orthogonal vector to v3 + ortho(1)=0.0d0 + ortho(2)=v3(3) + ortho(3)=-v3(2) + +!.. projected ONO angles in radians + pat_in(4)=get_ang(v2,v3,ortho) + pat_in(5)=get_ang(v1,v3,ortho) + + pat_in(6)=dabs(pat_in(5)-pat_in(4)) + +!.. account for rotational order of atoms + if (pat_in(4).le.pat_in(5)) then + pat_in(5)=2*pi-pat_in(4)-pat_in(6) + else + pat_in(4)=2*pi-pat_in(5)-pat_in(6) + endif + + pat_in(4)=rad2deg*pat_in(4)-offsets(2) + pat_in(5)=rad2deg*pat_in(5)-offsets(2) + pat_in(6)=rad2deg*pat_in(6)-offsets(2) + pat_in(7)=rad2deg*pat_in(7) + + call genANN_ctrans(pat_in) + + qint(:)=pat_in(:) + + contains + +!------------------------------------------------------------------- +! compute vector product n1 of vectors v1 x v2 + subroutine xprod(n1,v1,v2) + implicit none + + real*8 n1(3), v1(3), v2(3) + + n1(1) = v1(2)*v2(3) - v1(3)*v2(2) + n1(2) = v1(3)*v2(1) - v1(1)*v2(3) + n1(3) = v1(1)*v2(2) - v1(2)*v2(1) + + end subroutine + + +!------------------------------------------------------------------- +! compute scalar product of vectors v1 and v2: + real*8 function scalar(v1,v2,n) + implicit none + integer i, n + real*8 v1(*), v2(*) + + scalar=0.d0 + do i=1,n + scalar=scalar+v1(i)*v2(i) + enddo + + end function + +!------------------------------------------------------------------- +! compute norm of vector: + real*8 function norm(x,n) + implicit none + integer i, n + real*8 x(*) + + norm=0.d0 + do i=1,n + norm=norm+x(i)**2 + enddo + norm=sqrt(norm) + + end function + +!------------------------------------------------------------------- + + subroutine rot_trisec(tr,v1,v2,v3) + implicit none + + real*8 tr(3),v1(3),v2(3),v3(3) + + real*8 vrot(3) + real*8 rot_ax(3) + real*8 cos_phi,sin_phi + +! evaluate cos(-phi) and sin(-phi), where phi is the angle between +! tr and (1,0,0) + cos_phi=tr(1) + sin_phi=dsqrt(tr(2)**2+tr(3)**2) + + if (sin_phi.lt.1.0d-12) then + return + endif + +! determine rotational axis + rot_ax(1) = 0.0d0 + rot_ax(2) = tr(3) + rot_ax(3) = -tr(2) +! normalize + rot_ax=rot_ax/sin_phi + +! now the rotation can be done using Rodrigues' rotation formula +! v'=v*cos(p) + (k x v)sin(p) + k (k*v) (1-cos(p)) +! for v=tr k*v vanishes by construction: + +! check that the rotation does what it should + call rodrigues(vrot,tr,rot_ax,cos_phi,sin_phi) + if (dsqrt(vrot(2)**2+vrot(3)**2).gt.1.0d-12) then + write(6,*) "ERROR: BROKEN TRISECTOR" + stop + endif + tr=vrot + + call rodrigues(vrot,v1,rot_ax,cos_phi,sin_phi) + v1=vrot + call rodrigues(vrot,v2,rot_ax,cos_phi,sin_phi) + v2=vrot + call rodrigues(vrot,v3,rot_ax,cos_phi,sin_phi) + v3=vrot + + + end subroutine + +!------------------------------------------------------------------- + + subroutine rodrigues(vrot,v,axis,cos_phi,sin_phi) + implicit none + real*8 vrot(3),v(3),axis(3) + real*8 cos_phi,sin_phi + + real*8 ortho(3) + + call xprod(ortho,axis,v) + vrot = v*cos_phi + ortho*sin_phi + > + axis*scalar(axis,v,3)*(1-cos_phi) + + end subroutine + +!------------------------------------------------------------------- + + real*8 function get_ang(v,xaxis,yaxis) + implicit none +! get normalized [0:2pi) angle from vectors in the yz plane + real*8 v(3),xaxis(3),yaxis(3) + + real*8 phi + + real*8 pi + parameter (pi=3.141592653589793d0) + + phi=atan2(scalar(yaxis,v,3),scalar(xaxis,v,3)) + if (phi.lt.0.0d0) then + phi=2*pi+phi + endif + get_ang=phi + + end function + + end subroutine cart2int diff --git a/ctrans.f b/ctrans.f new file mode 100644 index 0000000..ca619dd --- /dev/null +++ b/ctrans.f @@ -0,0 +1,88 @@ +!------------------------------------------------------------------- +! Time-stamp: "2024-10-09 13:33:50 dwilliams" + + subroutine genANN_ctrans(pat_in) + implicit none + + include 'nnparams.incl' + include 'JTmod.incl' + + real*8 pat_in(maxnin) + + real*8 raw_in(maxnin),off_in(maxnin),ptrans_in(7) + real*8 r0 + real*8 a,b,xs,ys,xb,yb + + integer k + + off_in(1:7)=pat_in(1:7) + r0=offsets(1) + +! transform primitives +! recover raw distances from offset coords + do k=1,3 + raw_in(k)=off_in(k)+offsets(1) + enddo + + do k=1,3 + ptrans_in(k)=off_in(k) + enddo + +! rescale ONO angles + ptrans_in(4)=deg2rad*off_in(4) + ptrans_in(5)=deg2rad*off_in(5) + ptrans_in(6)=deg2rad*off_in(6) +! rescale umbrella + ptrans_in(7)=off_in(7)*deg2rad + +! compute symmetry coordinates + +! A (breathing) + a=(ptrans_in(1)+ptrans_in(2)+ptrans_in(3))/dsqrt(3.0d0) +! ES + call prim2emode(ptrans_in(1:3),xs,ys) +! EB + call prim2emode(ptrans_in(4:6),xb,yb) +! B (umbrella) + b=ptrans_in(7) + +! overwrite input with output + + pat_in(pat_index(1))=a ! 1 + pat_in(pat_index(2))=xs + pat_in(pat_index(3))=ys + pat_in(pat_index(4))=xb + pat_in(pat_index(5))=yb + pat_in(pat_index(6))=b +! totally symmetric monomials + pat_in(pat_index(7))=xs**2 + ys**2 ! 2 + pat_in(pat_index(8))=xb**2 + yb**2 ! 3 + pat_in(pat_index(9))=b**2 ! 9 + pat_in(pat_index(10))=xs*xb+ys*yb ! 4 +! S^3, B^3 + pat_in(pat_index(11))=xs*(xs**2-3*ys**2) ! 5 + pat_in(pat_index(12))=xb*(xb**2-3*yb**2) ! 6 +! S^2 B, S B^2 + pat_in(pat_index(13))=xb*(xs**2-ys**2) - 2*yb*xs*ys ! 7 + pat_in(pat_index(14))=xs*(xb**2-yb**2) - 2*ys*xb*yb ! 8 + + do k=11,14 + pat_in(pat_index(k))=tanh(0.1d0*pat_in(pat_index(k)))*10.0d0 + enddo + + contains + + subroutine prim2emode(prim,ex,ey) + implicit none +! Takes a 2D-vector prim and returns the degenerate modes x and y +! following our standard conventions. + + real*8 prim(3),ex,ey + + ex=(2.0d0*prim(1)-prim(2)-prim(3))/dsqrt(6.0d0) + ey=(prim(2)-prim(3))/dsqrt(2.0d0) + + end + + + end subroutine diff --git a/nnparams.incl b/nnparams.incl new file mode 100644 index 0000000..420cfc0 --- /dev/null +++ b/nnparams.incl @@ -0,0 +1,188 @@ +!**** Declarations + + real*8 pi,infty,zero + real*8 scan_res + real*8 hart2eV, eV2hart + real*8 hart2icm, icm2hart + real*8 eV2icm, icm2eV + real*8 deg2rad, rad2deg + integer maxneu,maxlay,maxtypes,maxtpar + integer maxpats + integer maxnin,maxnout,maxpout + integer maxwei,neucap,wbcap + integer maxxrmeta,xrcap + integer iinfty + integer iout,nnunit,perfunit,fitunit + integer ec_error,ec_read,ec_dim,ec_log + integer ec_dimrd + character*2 newline + character*8 stdfmt + character*8 nnldir + character*8 nntag + character*16 prim_tag + character*16 nnfdir,nnsdir + character*16 sline,asline,hline + character*16 mform,smform,miform + character*16 lrfmt,lifmt + character*32 nndmpfile,nnexpfile + character*32 nndatfile,nnreffile + character*32 sampfile,perfile + character*32 nnparfile,nnp10file + +!********************************************************** +!**** Parameters +!*** maxneu: max. number of neurons per hidden layer +!*** maxnin: max. number of neurons in input layer +!*** maxnout: max. number of neurons in output layer +!*** maxpout: max. number of values in output pattern +!*** maxlay: max. number of layers (always >2) +!*** maxtypes: max. number of neuron types +!*** maxtpar: max. number of parameters for each neuron type +!*** maxpats: max. number of learning patterns +!*** maxxrmeta: max. number of metadata-blocks in xranges + +!*** WARNING: maxnout should not be > maxneu, deriv-like structures +!*** assume so. + parameter (maxneu=150,maxnin=14,maxnout=15) + parameter (maxpout=15) + parameter (maxlay=3,maxtypes=2,maxtpar=1) + parameter (maxpats=50000) + + parameter (maxxrmeta=3) + +!********************************************************** +!**** Inferred Parameters +!*** maxwei: max. total number of weight matrix elements +!*** neucap: max. total number of neurons +!*** wbcap: max. total number of weights and biases +!*** xrcap: max. total number of used dimensions in xranges + + parameter (maxwei=(maxlay-3)*maxneu**2+maxneu*(maxnin+maxnout)) + parameter (neucap=(maxlay-2)*maxneu+maxnin+maxnout) + parameter (wbcap=maxwei+neucap) + parameter (xrcap=2+maxxrmeta) + +!*** WARNING: maxwei may fail for 2-layered networks +!*** if maxnin*maxnout is sufficiently large! + +!********************************************************** +!**** Numerical Parameters +!*** infty: largest possible double precision real value. +!*** iinfty: largest possible integer value. +!*** zero: sets what is considered an irrelevant difference +!*** in size. use for comarison of reals, to determine +!*** 'dangerously small' values, etc +!*** scan_res: maximum precision for geometric boundary algorithm + + ! 3.14159265358979323846264338327950... + parameter (pi=3.1415926536D0) + parameter (infty=huge(1.0D0),iinfty=huge(1)) + parameter (zero=1.0D-8,scan_res=1.0D-8) + +!********************************************************** +!**** Unit Conversion Parameters +!*** X2Y: convert from X to Y. +!*** +!**** !? currently inexact. FIX THIS. +!*** hart: hartree +!*** eV: electron volt +!*** icm: inverse centimeters (h*c/cm) +!**** +!*** deg: degree +!*** rad: radians + + parameter (hart2icm=219474.69d0) + parameter (hart2eV=27.211385d0) + parameter (eV2icm=hart2icm/hart2eV) + parameter (icm2hart=1.0d0/hart2icm) + parameter (eV2hart=1.0d0/hart2eV) + parameter (icm2eV=1.0d0/eV2icm) + parameter (deg2rad=pi/180.0d0) + parameter (rad2deg=1.0d0/deg2rad) + +!********************************************************** +!**** I/O Parameters +!*** iout: standard output for vranf error messages +!*** nnunit: temporary UNIT for misc. output files +!*** nnuit + [0..99] are reserved for futher +!*** unspecific misc. files. +!*** perfunit: UNIT for performance logfile +!*** fitunit: UNIT added to random positive integer +!*** identifying a single core fit UNIQUELY +!*** +!*** lrfmt: format for long real output +!*** lifmt: format for long integer output +!*** +!*** nndatfile: filename for DATA-files +!*** (without file extension) +!*** nnreffile: filename for reference DATA-blocks +!*** (without file extension) +!*** nnparfile: filename for best fitted parameters to be +!*** written on (without file extension) +!*** nnp10file: filename for the 10th percentile parameters to +!*** be written on (without file extension) +!*** nnexpfile: filename for modified neural network parameters +!*** (without file extension) +!*** sampfile: filename for displaying sampled points in +!*** configuration space +!*** nndmpfile: filename for dumping data point pairs +!*** perfile: filename for logged fitting performances. +!*** nntag: infix for various filenames to mark their origin +!*** program should end with a trailing '_' if nonempty. +!*** prim_tag: tag added to the '***' line of primitive par-files +!*** nnfdir: directory for dumping fit files +!*** nnsdir: directory for dumping scans. +!*** nnldir: directory for dumping logfiles for each fit + + parameter (nndatfile='DATA_ANN') + parameter (nnreffile='REF_ANN') + parameter (nnparfile='../nnfits/fit_pars') + parameter (nnp10file='../nnfits/fit_10p') + parameter (nnexpfile='../nnfits/exp_pars') + parameter (nndmpfile='../nnfits/fit_dump.dat') + parameter (sampfile='../scans/samples.dat') + parameter (perfile='../logs/performance.log') + parameter (nnfdir='../nnfits/',nnsdir='../scans/') + parameter (nnldir='../logs/') + parameter (nntag='',prim_tag=' Time-stamp: " "') + + parameter (lrfmt='(ES20.12)',lifmt='(I)') + + parameter (iout=6,perfunit=700,nnunit=800,fitunit=8000) + +!********************************************************** +!**** Debugging Parameters +!*** sline: separation line +!*** asline: alternative sep. line +!*** hline: simple horizontal line +!*** newline: a single blank line +!*** mform: standard form for matrix output +!*** miform: standard form for integer matrix output +!*** smform: shortened form for matrix output +!*** stdfmt: standard format for strings + + parameter (sline='(75("*"))',asline='(75("#"))') + parameter (hline='(75("-"))') + parameter (newline='()') + parameter (mform='(5ES12.4)',smform='(5ES10.2)') + parameter (miform='(5I12)') + parameter (stdfmt='(A)') + +!********************************************************** +!**** Error Codes +!*** Codes should be powers of 2. Binary representation of return value +!*** should correspond to all exceptions invoked. ec_error should never +!*** be invoked with any other. +!*** +!*** ec_error: generic error (catch-all, avoid!) +!*** ec_read: parsing error during les() +!*** ec_dim: dimensioning error +!*** ec_log: logic error +!*** +!**** Inferred error codes +!*** ec_dimrd: ec_dim+ec_read + + + parameter (ec_error=1,ec_read=2,ec_dim=4,ec_log=8) + + parameter (ec_dimrd=ec_dim+ec_read) diff --git a/tab2genetic.f b/tab2genetic.f new file mode 100644 index 0000000..fa28fa5 --- /dev/null +++ b/tab2genetic.f @@ -0,0 +1,48 @@ + program tab2genetic + use iso_fortran_env, only: error_unit + implicit none + include 'nnparams.incl' + include 'JTmod.incl' + + integer, parameter :: infile=100 + character(len=2048) line + + real*8 cart(3,4),qint(maxnin) +! useless crap written in table1 + real*8 tmp + + character(len=1024) binary,infile_name + character(len=3) lfmt,ofmt + parameter (lfmt='(A)',ofmt='(A)') + integer, parameter :: uerr=error_unit, uout=6 + + integer j + + call getarg(0,binary) + + if (iargc().lt.1) then + write(uerr,ofmt) 'ERROR: Missing first argument FILENAME.' + write(uerr,ofmt) trim(binary) + > // ' FILENAME' + stop + endif + call getarg(1,infile_name) + write(uerr,ofmt) ' Input file: '//trim(infile_name) + open(infile,file=trim(infile_name),status='old',action='read') + + do + read(infile,lfmt,err=403,end=404) line + cart(:,1)=0 + read(line,*,err=405,end=403) tmp, cart(1:3,2:4) + call cart2int(cart,qint) +! order: xs,ys,xb,yb,a,b + write(uout,'(6F20.6)') (qint(pat_index(j)), j=2,5), + > qint(pat_index(1)), qint(pat_index(6)) + 403 cycle + 404 exit + 405 write(uerr,*) 'WARNING: MALFORMED LINE: "' + > //trim(line)//'"' + enddo + + close(infile) + end program