From ccf300a77e7d188a4783eb70abbd380ba62efb67 Mon Sep 17 00:00:00 2001 From: jean paul nshuti Date: Thu, 9 Oct 2025 10:05:57 +0200 Subject: [PATCH] Add the model for planar geometry of NH3+ --- Makefile | 177 ++++++++ src/model/genetic_param_no3.f90 | 1 + src/model/invariants_no3.f90 | 82 ++++ src/model/model_no3.f90 | 469 +++++++++++++++++++++ src/model/nnadia_no3.f | 63 +++ src/model/nncoords_no3.f90 | 716 ++++++++++++++++++++++++++++++++ src/model/red_invariants.f90 | 72 ++++ 7 files changed, 1580 insertions(+) create mode 100644 Makefile create mode 120000 src/model/genetic_param_no3.f90 create mode 100644 src/model/invariants_no3.f90 create mode 100644 src/model/model_no3.f90 create mode 100644 src/model/nnadia_no3.f create mode 100644 src/model/nncoords_no3.f90 create mode 100644 src/model/red_invariants.f90 diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..96e973b --- /dev/null +++ b/Makefile @@ -0,0 +1,177 @@ +###################################################################### +# GNU MAKEFILE # +###################################################################### +# There may be some trickery involved in this file; most of which is # +# hopefully explained sufficiently. This makefile is unlikely to # +# work with non-GNU make utilities. # +# # +###################################################################### +SHELL = /bin/bash +# clear out suffix list +.SUFFIXES : +# declare suffixes allowed to be subject to implicit rules. +.SUFFIXES : .f .o .f90 + +current_version=4.0c +tarname=NH3-Plus + + +src = ./src/ +build = ./obj/ +bin = ./bin/ + +ldir = $(src)lib/ +pdir = $(src)parser/ + + +# DEFINE THE PATH + +VPATH = $(src) +VPATH += $(ldir) +VPATH += $(pdir) + + +# extra dir +extra_dirs = logs nnfits scans +gincl = -I$(src) -I$(pdir) -I$(ldir) + +# GFORTRAN COMPILER +FC := gfortran +GFFLAGS = -O2 -fopenmp -fbackslash -fmax-errors=5 -Wall -Werror +GFFLAGS += -std=legacy -cpp -J$(build) +GFFLAGS += $(gincl) +#GFFLAGS += -Wall -Wextra -Werror -Wno-error=integer-division -Wno-error=conversion \ + -Wno-error=intrinsic-shadow +GLLAPCK = -llapack + +DBGFFLAGS = -O0 -fcheck=bounds -fcheck=do -fcheck=mem -fcheck=pointer -p -debug all + +#cGINCL = -I -I$(pdir) -I$(ldir) + +GFFLAGS += $(GINCL) +GFFLAGS += $(GLLAPCK) +DBGFFLAGS += $(GLLAPCK) + +#IFORT COMPILER + +IFC := ifort +IFFLAGS = -O3 -qopenmp -qmkl -fpp -stand f95 -warn all -warn noerrors \ + -fp-model precise -traceback -assume byterecl +IFFLAGS += -g -diag-disable=10448 +DBGIFC = -O0 -g -check all -traceback -fpp -warn all -fp-stack-check + +LDFLAGS = -llapack + +# CHOSE THE COMPILER +########################### +COMPILER = $(FC) + +ifeq ($(COMPILER),IFC) + FC = $(IFC) + FFLAGS = $(IFFLAGS) + DBGFLAGS = $(DBGIFC) +else + #FC = $(FC) + FFLAGS = $(GFFLAGS) + DBGFLAGS = $(DBGFFLAGS) +endif +###################################################################### +### Functions +###################################################################### + +# return 'true' if file $(1) exists, else return '' +file_exists = $(shell if [ -f $(1) ]; then echo "true"; fi) +# return 'true' if file $(1) doesn't exist, else return '' +file_lost = $(shell if [ ! -f $(1) ]; then echo "true"; fi) + +###################################################################### +### Objects +###################################################################### + +# Objects solely relying on ANN parameters +ann_objects = backprop.o ff_neunet.o neuron_types.o \ + axel.o nnmqo.o scans.o nnread.o + +# Objects relying on both genetic & ANN params +genann_objects = geNNetic.o iNNterface.o puNNch.o \ + mkNN.o error.o long_io.o parse_errors.o parser.o + +# Objects depending on model +#mod_objects = nnmodel.o nncoords.o +pln_objects = invariants_no3.o nncoords_no3.o genetic_param_no3.o model_no3.o nnadia_no3.o +#pyr_objects = invariants_nh3.o nncoords_nh3.o fit_genetic_Umb.o model_dipole_nh3.o dip_nh3plus_model.o +#mod_objects = invariants_nh3.o nncoords_nh3.o fit_genetic.o \ + new_model_dipole.o dip_nh3plus_model.o + +# Objects belonging to internal library +lib_objects = misc.o ran_gv.o choldc.o diag.o \ + qsort.o dmatrix.o imatrix.o strings.o ran.o \ + fileread.o keyread.o long_keyread.o + +# Objects from hell; if it exists, don't compile it +hell_objects = $(ldir)random.o +hell_flags = -O3 +hell_flags += -llapack + +# All objects (one would want to update) +objects = $(addprefix $(build), $(ann_objects) $(genann_objects) $(pln_objects) \ + $(lib_objects)) +#umb_objects = $(addprefix $(build), $(ann_objects) $(genann_objects) $(pyr_objects) \ + $(lib_objects)) +#surface_obj = $(addprefix $(build), invariants_nh3.o nncoords_nh3.o nnread.o \ + try_param.o model_dipole_nh3.o) +# Name of the main program +main = genANN + +###################################################################### +### Include files +###################################################################### + +ann_include = nncommon.incl nndbg.incl nnparams.incl +gen_include = common.incl params.incl keylist.incl errcat.incl +mod_include = JTmod.incl only_model.incl dip_planar_genetic.incl + +# Misc. files generated during compilation +trash = *__genmod* *~ *\# *.g .mod + +$(main) : $(objects) +ifeq ($(call file_exists,$(hell_objects)),) + $(MAKE) hell +endif + $(FC) $(FFLAGS) $(hell_objects) $(objects) $(LDFLAGS) -o $(bin)$(main) +#$(FFLAGS) + +#pyram_ANN : $(umb_objects) +# $(FC) $(FFLAGS) $(hell_objects) $(umb_objects) $(LDFLAGS) -o $(bin)pyram_ANN +# compile a library for the Diabatic Dipole + +#libdipole_surface.a : $(surface_obj) $(build)Diabatic_Dipole.o +# ar rcs $(bin)libdipole_surface.a $(surface_obj) $(build)Diabatic_Dipole.o +# EXPLICIT RULE FOR COMPILING + +$(build)%.o : %.f + $(FC) -c $(FFLAGS) $< -o $@ +$(build)%.o : %.f90 + $(FC) -c $(FFLAGS) $< -o $@ + +# add dependencies to include files +$(ann_objects) : $(ann_include) $(lib_objects) +$(genann_objects) : $(ann_include) $(gen_include) $(lib_objects) +$(mod_objects) : $(ann_include) $(gen_include) $(mod_include) $(lib_objects) +# overrule standard compilation method for hellfiles. +$(hell_objects) : override FFLAGS=$(hell_flags) + +nnmqo.o : axel.o + +# making + +.PHONY : clean dirs neat hell pyram_ANN all +clean: + $(RM) $(objects) $(hell_objects) $(trash) $(bin)/$(main) +dirs: + @mkdir -p $(build) $(bin) $(extra_dirs) +neat: + $(RM) -f $(TRASH) +hell : $(hell_objects) + +all : clean hell $(main) pyram_ANN libdipole_surface.a diff --git a/src/model/genetic_param_no3.f90 b/src/model/genetic_param_no3.f90 new file mode 120000 index 0000000..50a5cb2 --- /dev/null +++ b/src/model/genetic_param_no3.f90 @@ -0,0 +1 @@ +../../../../Genetic/NO3/Dipole_NO3/Fit_stretch_Latest/fit_genric_bend_no3.f90 \ No newline at end of file diff --git a/src/model/invariants_no3.f90 b/src/model/invariants_no3.f90 new file mode 100644 index 0000000..d78f08b --- /dev/null +++ b/src/model/invariants_no3.f90 @@ -0,0 +1,82 @@ +Module invariants_mod + implicit none + contains + !---------------------------------------------------- + subroutine invariants(a,xs,ys,xb,yb,b,inv) + implicit none + !include "nnparams.incl" + double precision, intent(in) :: a, xs, ys, xb, yb, b + double precision, intent(out) :: inv(3) + double precision:: invar(24) + complex(8) :: q1, q2 + LOGICAL,PARAMETER:: debg =.FALSE. + integer :: i + ! express the coordinate in complex + + q1 = dcmplx(xs, ys) + q2 = dcmplx(xb, yb) + + ! compute the invariants + invar(24) = a + invar(23) =b**2 + + ! INVARIANTS OF KIND II + !------------------------ + + invar(1) = dreal( q1 * conjg(q1) ) ! r11 + invar(2) = dreal( q1 * conjg(q2) ) ! r12 + invar(3) = dreal( q2 * conjg(q2) ) ! r22 + invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2 + + + !INVATIANTS OF KIND III + !------------------------ + + invar(5) = dreal( q1 * q1 * q1 ) ! r111 + invar(6) = dreal( q1 * q1 * q2 ) ! r112 + invar(7) = dreal( q1 * q2 * q2 ) ! r122 + invar(8) = dreal( q2 * q2 * q2 ) ! r222 + invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2 + invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2 + invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2 + invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222 + + ! INVARIANTS OF KIND IV + !------------------------- + + invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 )) + invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 )) + invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 )) + invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 )) + + ! INVARIANTS OF KIND V + !---------------------- + + invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 )) + invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 )) + invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 )) + invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 )) + invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 )) + invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 )) + + ! the only non zero invariant for bend pure cuts + + inv(1) = invar(3) + inv(2) = invar(8) + inv(3) = invar(12) + + if (debg) then + write(*,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) + write(*,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) + write(*,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) + write(*,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) + + write(*,*)"THE INPUT COORDINATE IN COMPLEX REPRES" + write(*,*)"---------------------------------------" + write(*,*)"xs =",dreal(q1), "ys=",dimag(q1) + endif + ! modify the invariants to only consider few of them + ! + !invar(13:22)=0.0d0 + end subroutine invariants +end module invariants_mod diff --git a/src/model/model_no3.f90 b/src/model/model_no3.f90 new file mode 100644 index 0000000..8fadd6e --- /dev/null +++ b/src/model/model_no3.f90 @@ -0,0 +1,469 @@ +module diabmodel + !use dim_parameter,only:qn,ndiab,pst + use iso_fortran_env, only: idp => int32, dp => real64 + !use accuracy_constants, only:dp,idp + use dip_param, only: init_dip_planar_data, p,pst,np + + implicit none + include "nnparams.incl" + integer(idp),parameter:: ndiab=5 + logical :: debug=.false. + contains + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + subroutine diab_x(q,e,nn_out) + real(dp),intent(in)::q(maxnin) + real(dp),intent(out)::e(ndiab,ndiab) + real(dp),intent(inout):: nn_out(maxnout) + integer(idp) id,i,j, ii, non_zer_p + real(dp) xs,xb,ys,yb,a,b,ss,sb,v3(8),v2(6) + xs=q(5) + ys=q(6) + xb=q(7) + yb=q(8) + a=q(4) + b=q(9) + + call init_dip_planar_data() + non_zer_p = count(p /= 0.0d0) + + ii=1 + do i=1,np-2 + if (p(i) .ne. 0) then + p(i) =p(i)*(1.0_dp + 1.0d-2 + nn_out(ii)) + ii=ii+1 + else + p(i)=p(i) + endif + enddo + ss=xs**2+ys**2 ! totaly symmetric term + sb=xb**2+yb**2 + v2(1)=xs**2-ys**2 + v2(2)=xb**2-yb**2 + v2(3)=xs*xb-ys*yb + v2(4)=2*xs*ys + v2(5)=2*xb*yb + v2(6)=xs*yb+xb*ys + + v3( 1) = xs*(xs**2-3*ys**2) + v3( 2) = xb*(xb**2-3*yb**2) + v3( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys + v3( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb + v3( 5) = ys*(3*xs**2-ys**2) + v3( 6) = yb*(3*xb**2-yb**2) + v3( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys + v3( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb + + + + e=0.0d0 + id=1 ! 1 + e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 param + id=id+1 ! 2 + e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 p + e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + id =id+1 ! 3 + e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 p + e(5,5)=e(5,5)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb + + + id=id+1 ! 4 + ! order 2 + e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 5 p + +p(pst(1,id)+2)*(xs*xb-ys*yb) + id =id+1 ! 5 + e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + id =id+1 ! 6 + e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + e(5,5)=e(5,5)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + ! order 3 + id=id+1 ! 7 + e(1,1)=e(1,1)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb *sb ! 2 param + id=id+1 ! 8 + e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb *sb ! 2 p + e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + id =id+1 ! 9 + e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb ! 2 p + e(5,5)=e(5,5)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + + + ! W and Z term of E1 + ! order 0 + id=id+1 ! 10 + e(2,2)=e(2,2)+p(pst(1,id)) + e(3,3)=e(3,3)-p(pst(1,id)) + !e(2,3)=e(2,3) + + ! order 1 + id=id+1 ! 11 ! 2 param + e(2,2)=e(2,2)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb + e(3,3)=e(3,3)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb) + e(2,3)=e(2,3)- p(pst(1,id))*ys -p(pst(1,id)+1)*yb + ! order 2 + id=id+1 ! 12 ! 3p + do i=1,3 + e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i) + e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i) + e(2,3)=e(2,3)+ p(pst(1,id)+(i-1))*v2(i+3) + enddo + ! order 3 + id=id+1 ! 13 ! 8 param + do i=1,4 + e(2,2)=e(2,2)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + e(3,3)=e(3,3)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + e(2,3)=e(2,3)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i+4) + enddo + + ! try the testing of higher order terms + !e(2,3)=e(2,3)- p(pst(1,id))*ys*ss +p(pst(1,id)+1)*ss*2*xs*ys + + + + + ! W and Z for E2 + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + id=id+1 ! 14 + e(4,4)=e(4,4)+p(pst(1,id)) + e(5,5)=e(5,5)-p(pst(1,id)) + e(4,5)=e(4,5) + + ! order 1 + id=id+1 ! 2 param 15 + e(4,4)=e(4,4)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb + e(5,5)=e(5,5)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb) + e(4,5)=e(4,5)- p(pst(1,id))*ys-p(pst(1,id)+1)*yb + ! order 2 + id=id+1 ! 16 ! 3p + do i=1,3 + e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i) + e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i) + e(4,5)=e(4,5)+ p(pst(1,id)+(i-1))*v2(i+3) + enddo + + ! order 3 + id=id+1 ! 17 ! 8 param + do i=1,4 + e(4,4)=e(4,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + e(5,5)=e(5,5)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + e(4,5)=e(4,5)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i+4) + enddo + + !e(4,4) = e(4,4)+ p(pst(1,id))*xs*ss + p(pst(1,id)+1)*xb*sb + !e(5,5)=e(5,5)- (p(pst(1,id))*xs*ss + p(pst(1,id)+1)*xb*sb) + !e(4,5)=e(4,5)- p(pst(1,id))*ys*ss -p(pst(1,id)+1)*yb*sb + + ! E1 X E2 + ! WW and ZZ + + + id =id+1 ! 18 + e(2,4)=e(2,4)+p(pst(1,id))*b + e(3,5)=e(3,5)-p(pst(1,id))*b + + ! ORDER 1 + id=id+1 ! 19 ! 6 parama + e(2,4)=e(2,4)+b*((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*xb) + e(3,5)=e(3,5)+b*((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb) + e(2,5)=e(2,5)+b*((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*yb) + e(3,4)=e(3,4)+b*((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb) + ! order 2 + id=id+1 ! 20 + + do i=1,3 ! 9 param + e(2,4)=e(2,4)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i) + e(3,5)=e(3,5)+b*(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i) + e(2,5)=e(2,5)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3) + e(3,4)=e(3,4)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i+3) + enddo + + ! pseudo A2 & E1 + ! ################################################## + !################################################### + ! order 0 + id=id+1 ! 1 param ! 21 + + e(1,3)=e(1,3)+b*(p(pst(1,id))) + + ! order 1 + id = id +1 ! 22 + e(1,2)=e(1,2)-b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) + e(1,3)=e(1,3)+b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) + + ! order 2 + id=id+1 ! 23 + + e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)) + e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb)) + + ! COUPLING OF A2 WITH E2 + + !########################################################################################################## + + ! order 0 + id =id+1 !24 + e(1,5)=e(1,5)+p(pst(1,id)) + + ! order 1 + id = id +1 ! 25 + e(1,4)=e(1,4)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) + e(1,5)=e(1,5)+(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) + + ! order 2 + id=id+1 ! 26 + + e(1,4)=e(1,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys) + e(1,5)=e(1,5)+p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb) + + ! order 3 + id=id+1 ! 27 ! 8 param + do i=1,4 + e(1,4)=e(1,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4) + e(1,5)=e(1,5)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + + enddo + + + call copy_2_lower_triangle(e) + !temp = e(2,2) + !e(2,2)=e(3,3) + !e(3,3)=temp + if (debug) then + do i=1,ndiab + write(34,'(5f14.6)') (e(i,j),j=1,ndiab) + enddo + write(34,*)"" + endif + + end subroutine diab_x + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! THE Y COMPONENT OF DIPOLE + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine diab_y(q,e,nn_out) + real(dp),intent(in)::q(maxnin) + real(dp),intent(out)::e(ndiab,ndiab) + real(dp),intent(inout):: nn_out(maxnout) + integer(idp) id,i,j, ii, non_zer_p + real(dp) xs,xb,ys,yb,a,b,ss,sb,v3(8),v2(6) + xs=q(5) + ys=q(6) + xb=q(7) + yb=q(8) + a=q(4) + b=q(9) + + call init_dip_planar_data() + non_zer_p = count(p /= 0.0d0) + + ii=1 + do i=1,np-2 + if (p(i) .ne. 0) then + p(i) =p(i)*(1.0_dp + 1.0d-2 + nn_out(ii)) + ii=ii+1 + else + p(i)=p(i) + endif + enddo + ss=xs**2+ys**2 ! totaly symmetric term + sb=xb**2+yb**2 + v2(1)=xs**2-ys**2 + v2(2)=xb**2-yb**2 + v2(3)=xs*xb-ys*yb + v2(4)=2*xs*ys + v2(5)=2*xb*yb + v2(6)=xs*yb+xb*ys + + + v3( 1) = xs*(xs**2-3*ys**2) + v3( 2) = xb*(xb**2-3*yb**2) + v3( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys + v3( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb + v3( 5) = ys*(3*xs**2-ys**2) + v3( 6) = yb*(3*xb**2-yb**2) + v3( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys + v3( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb + + e=0.0d0 + ! V-term + + id=1 ! 1 + ! order 1 + e(1,1)=e(1,1)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb + id=id+1 ! 2 + e(2,2)=e(2,2)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb + e(3,3)=e(3,3)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb + id =id+1 ! 3 + e(4,4)=e(4,4)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb + e(5,5)=e(5,5)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb + + id=id+1 ! 4b*( + e(1,1)=e(1,1)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)) + id =id+1 ! 5 + e(2,2)=e(2,2)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)) + e(3,3)=e(3,3)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)) + id=id+1 ! 6 + e(4,4)=e(4,4)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)) + e(5,5)=e(5,5)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)) + ! order 3 + id=id+1 ! 7 + + e(1,1)=e(1,1)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb + id=id+1 ! 2 + e(2,2)=e(2,2)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb + e(3,3)=e(3,3)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb + id =id+1 ! 3 + e(4,4)=e(4,4)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb + e(5,5)=e(5,5)+p(pst(1,id))*ys*ss + p(pst(1,id)+1)*yb*sb + + ! W and Z of E1 + ! order 0 + id=id+1 ! 10 + e(2,3)=e(2,3)+p(pst(1,id)) + ! order 1 + id=id+1 ! + e(2,2)=e(2,2)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb + e(3,3)=e(3,3)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb + e(2,3)=e(2,3)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb + ! order 2 + id=id+1 ! 12 + do i=1,3 + e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i+3) + e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i+3) + e(2,3)=e(2,3)-p(pst(1,id)+(i-1))*v2(i) + enddo + + id=id+1 ! 8 + do i=1,4 + e(2,2)=e(2,2)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4) + e(3,3)=e(3,3)-(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4) + e(2,3)=e(2,3)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + enddo + + + !! W and Z of E2 + !!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! order 0 + id=id+1 ! 14 + e(4,5)=e(4,5)+p(pst(1,id)) + ! order 1 + id=id+1 ! 15 + e(4,4)=e(4,4)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb + e(5,5)=e(5,5)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb + e(4,5)=e(4,5)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb + ! order 2 + id=id+1 ! 16 + do i=1,3 + e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i+3) + e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i+3) + e(4,5)=e(4,5)-p(pst(1,id)+(i-1))*v2(i) + enddo + + id=id+1 ! 17 + do i=1,4 + e(4,4)=e(4,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4) + e(5,5)=e(5,5)-(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3(i+4) + e(4,5)=e(4,5)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3(i) + enddo + ! PSEUDO JAHN-TELLER E1 AND E2 + + !ORDER 0 + id=id+1 ! 18 + + e(2,5)=e(2,5)+p(pst(1,id))*b + e(3,4)=e(3,4)+p(pst(1,id))*b + ! order 1 + + id=id+1 + e(2,4)=e(2,4)+b*((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb) + e(3,5)=e(3,5)+b*((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*yb) + e(2,5)=e(2,5)+b*((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb) + e(3,4)=e(3,4)+b*((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(+p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*xb) + + ! order 2 + id=id+1 + + e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)*b + e(3,5)=e(3,5)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)*b + e(2,5)=e(2,5)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i)*b + e(3,4)=e(3,4)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)*b + + ! no order 3 + !!!!!!!!!!!!!!!! + + ! the coupling A2 & E1 + ! ##################### + ! order 0 + + id=id+1 + e(1,2)=e(1,2)+b*(p(pst(1,id))) + ! order 1 + + id=id+1 + e(1,2)=e(1,2)-b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) + e(1,3)=e(1,3)-b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) + + ! order 2 + id=id+1 + e(1,2)=e(1,2)-b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb)) + e(1,3)=e(1,3)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ & + +p(pst(1,id)+2)*(xs*yb+xb*ys)) + + ! COUPLING OF A2 WITH E2 + !####################################################################################### + !############################################################################### + ! order 0 + + id = id+1 + e(1,4)=e(1,4)+p(pst(1,id)) + ! order 1 + + id=id+1 + e(1,4)=e(1,4)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb) + e(1,5)=e(1,5)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) + + ! order 2 + id=id+1 + e(1,4)=e(1,4)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & + +p(pst(1,id)+2)*(xs*xb-ys*yb)) + e(1,5)=e(1,5)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ & + p(pst(1,id)+2)*(xs*yb+xb*ys)) + + !write(*,*)'idy=',id + call copy_2_lower_triangle(e) + + if (debug) then + do i=1,ndiab + write(*,'(5f14.6)') (e(i,j),j=1,ndiab) + enddo + write(*,*)"" + endif +end subroutine diab_y + + + subroutine copy_2_lower_triangle(mat) + real(dp), intent(inout) :: mat(:, :) + integer :: m, n +! write lower triangle of matrix symmetrical + do n=1,size(mat,1) + do m=n,size(mat,1) + mat(m,n)=mat(n,m) + enddo + enddo + end subroutine copy_2_lower_triangle + +end module diabmodel diff --git a/src/model/nnadia_no3.f b/src/model/nnadia_no3.f new file mode 100644 index 0000000..71c1fc1 --- /dev/null +++ b/src/model/nnadia_no3.f @@ -0,0 +1,63 @@ + subroutine nninit_mod() + implicit none + + include 'nnparams.incl' + end subroutine + + + subroutine nnadia(coords,coeffs,adiaoutp) + USE diabmodel, only: diab_x + implicit none + ! returns THE DIABATIC MATRIX UPPER MATRIX OF BOTH X AND Y COMPONENT OF DIPOLE +! +! coords: vector containing symmetry-adapted coordinates. +! coeffs: vector containing coeffs of diab dipole matrix + +! adiaoutp: upper diab matrix of dipole +! x component come first and then y component +! +! dmat_x & dmat_y: analytic model of dipole, x and y +! +! matdim: dimension of hamiltoniam matrix +! mkeigen: integer: no eigenvectors if =0 +! eigenvs: dummy storage for eigenvectors + + include 'nnparams.incl' + include 'JTmod.incl' + integer ndiab + parameter ndiab=5 + DOUBLE PRECISION coords(maxnin),coeffs(maxnout) + DOUBLE PRECISION adiaoutp(maxpout),dmat_x(ndiab,ndiab) + !DOUBLE PRECISION dmat_y(ndiab,ndiab) + !integer i,j,ii + call diab_x(coords,dmat_x,coeffs) + !call diab_y(coords,dmat_y,coeffs) + + ! rewrite the diabatic matrix into 1D array of adiaoutp + !ii=1 + !do i=1,ndiab + !do j=i,ndiab + ! adiaoutp(ii)=dmat_x(i,j) + ! adiaoutp(ii+10) = -1*dmat_y(i,j) + ! ii=ii+1 + !enddo + !enddo + adiaoutp(1) = dmat_x(1,1) + adiaoutp(2) = dmat_x(2,1) + adiaoutp(3) = dmat_x(2,2) + adiaoutp(4) = dmat_x(3,1) + adiaoutp(5) = dmat_x(3,2) + adiaoutp(6) = dmat_x(3,3) + adiaoutp(7) = dmat_x(4,1) + adiaoutp(8) = dmat_x(4,2) + adiaoutp(9) = dmat_x(4,3) + adiaoutp(10) = dmat_x(4,4) + adiaoutp(11) = dmat_x(5,1) + adiaoutp(12) = dmat_x(5,2) + adiaoutp(13) = dmat_x(5,3) + adiaoutp(14) = dmat_x(5,4) + adiaoutp(15) = dmat_x(5,5) + + !write(*,*) dmat_x(1,1) + END SUBROUTINE + diff --git a/src/model/nncoords_no3.f90 b/src/model/nncoords_no3.f90 new file mode 100644 index 0000000..ab0a9ad --- /dev/null +++ b/src/model/nncoords_no3.f90 @@ -0,0 +1,716 @@ + +!-------------------------------------------------------------------------------------- + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % SUBROUTINE CTRANS(...) +! % +! % M. Vossel 21.03.2023 +! % +! % Routine to transform symmetryinput coordinates to symmetrized +! % coordinates. Distances Are discribet by Morse coordinates or +! % TMC depending on Set Parameters in the Genetic Input. +! % +! % input variables +! % q: +! % q(1): H1x +! % q(2): y +! % q(3): z +! % q(4): H2x +! % q(5): y +! % q(6): z +! % q(7): H3x +! % q(8): y +! % q(9): z +! +! +! +! % Internal variables: +! % t: primitive coordinates (double[qn]) +! % t(1): +! % t(2): +! % t(3): +! % t(4): +! % t(5): +! % t(6): +! % t(7): +! % t(8): +! % t(9): +! % t: dummy (double[qn]) +! % p: parameter vector +! % npar: length of parameter vector +! % +! % Output variables +! % s: symmetrized coordinates (double[qn]) +! % s(1): CH-symetric streatch +! % s(2): CH-asymetric streatch-ex +! % s(3): CH-asymetric streatch-ey +! % s(4): CH-bend-ex +! % s(5): CH-bend-ey +! % s(6): CH-umbrella +! % s(7): CH-umbrella**2 +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +module ctrans_mod + implicit none + include 'only_model.incl' + include 'nnparams.incl' +! precalculate pi, 2*pi and angle to radian conversion + double precision, parameter :: pii = 4.00d0*datan(1.00d0) + double precision, parameter :: pi2 = 2.00d0*pii + double precision, parameter :: ang2rad = pii/180.00d0 +! precalculate roots + double precision, parameter:: sq2 = 1.00d0/dsqrt(2.00d0) + double precision, parameter:: sq3 = 1.00d0/dsqrt(3.00d0) + double precision, parameter:: sq6 = 1.00d0/dsqrt(6.00d0) +! change distances for equilibrium + double precision, parameter :: dchequi = 2.344419d0 + +contains + subroutine ctrans(q) + !use dim_parameter, only: qn + use invariants_mod, only: invariants + integer k !running indices + double precision, intent(inout) :: q(maxnin) !given coordinates + double precision :: s(maxnin) !output coordinates symmetry adapted and scaled + double precision :: t(maxnin) !output coordinates symmetry adapted but not scaled + ! ANN Variables + !double precision, optional, intent(out) :: invariants(:) + ! kartesian coordianates copy from MeF+ so substitute c by n and removed f + double precision ch1(3), ch2(3), ch3(3), c_atom(3) + double precision nh1(3), nh2(3), nh3(3) + double precision zaxis(3), xaxis(3), yaxis(3) + double precision ph1(3), ph2(3), ph3(3) + ! primitive coordinates + double precision dch1, dch2, dch3 !nh-distances + double precision umb !Umbrella Angle from xy-plane + + ! Symmetry coordinates + double precision aR !a1-modes H-Dist., + double precision exR, exAng !ex components H-Dist., H-Ang. + double precision eyR, eyAng !ey components H-Dist., H-Ang. + double precision inv(3) + ! debugging + logical, parameter :: dbg = .false. + + ! initialize coordinate vectors + s = 0.0d0 + t = 0.0d0 + + ! write kartesian coords for readability + c_atom(1:3) = q(1:3) ! N-atom at origin + do k = 1, 3 + ch1(k) = q(k + 3) + ch2(k) = q(k + 6) + ch3(k) = q(k + 9) + end do + q=0.d0 + + ! construct z-axis + nh1 = normalized(ch1) + nh2 = normalized(ch2) + nh3 = normalized(ch3) + zaxis = create_plane(nh1, nh2, nh3) + + ! calculate bonding distance + dch1 = norm(ch1) + dch2 = norm(ch2) + dch3 = norm(ch3) + + ! construct symmertic and antisymmetric strech + aR = symmetrize(dch1 - dchequi, dch2 - dchequi, dch3 - dchequi, 'a') + exR = symmetrize(dch1, dch2, dch3, 'x') + eyR = symmetrize(dch1, dch2, dch3, 'y') + + ! construc x-axis and y axis + ph1 = normalized(project_point_into_plane(nh1, zaxis, c_atom)) + xaxis = normalized(ph1) + yaxis = xproduct(zaxis, xaxis) ! right hand side koordinates + + ! project H atoms into C plane + ph2 = normalized(project_point_into_plane(nh2, zaxis, c_atom)) + ph3 = normalized(project_point_into_plane(nh3, zaxis, c_atom)) + + call construct_HBend(exAng, eyAng, ph1, ph2, ph3, xaxis, yaxis) + umb = construct_umbrella(nh1, nh2, nh3, zaxis) + + ! set symmetry coordinates and even powers of umbrella + !s(1) = dch1-dchequi!aR + !s(2) = dch2-dchequi!exR + !s(3) = dch3-dchequi!eyR + ! call invariants and get them + ! 24 invariants + + call invariants(0.0d0,exR,eyR,exAng,eyAng,umb,inv) + q(1:3)=inv(1:3) + q(4) = aR + q(5) = exR + q(6) = eyR + q(7) = exAng + q(8) = -1.0d0*eyAng + q(9) = umb + ! pairwise distances as second coordinate set + !call pair_distance(q, t(1:6)) + + !if (dbg) write (6, '("sym coords s=",9f16.8)') s(1:qn) + !if (dbg) write (6, '("sym coords t=",9f16.8)') t(1:qn) + !if (present(invariants)) then + ! call get_invariants(s, invariants) + !end if + ! RETURN Q AS INTERNAL COORD + end subroutine ctrans +! subroutine ctrans(q) +! use invariants_mod, only: invariants +! implicit none +! double precision, intent(inout):: q(maxnin) +! double precision:: invar(24) +! double precision:: a,b,esx,esy,ebx,eby +! a=q(1) +! esx=q(2) +! esy=q(3) +! ebx=q(4) +! eby=q(5) +! b=q(6) +! call invariants(a,esx,esy,ebx,eby,b,invar) +! +! q(1:24)=invar(1:24) +! q(25)=esx +! q(26)=esy +! q(27)=ebx +! q(28)=eby +! q(29)=b +! end subroutine ctrans + subroutine pair_distance(q, r) + double precision, intent(in) :: q(9) + double precision, intent(out) :: r(6) + double precision :: atom(3, 4) + integer :: n, k, count + + !atom order: H1 H2 H3 N + atom(:, 1:3) = reshape(q, [3, 3]) + atom(:, 4) = (/0.00d0, 0.00d0, 0.00d0/) + + ! disntace order 12 13 14 23 24 34 + count = 0 + do n = 1, size(atom, 2) + do k = n + 1, size(atom, 2) + count = count + 1 + r(count) = sqrt(sum((atom(:, k) - atom(:, n))**2)) + end do + end do + end subroutine pair_distance + + function morse_and_symmetrize(x,p,pst) result(s) + double precision, intent(in),dimension(3) :: x + double precision, intent(in),dimension(11) :: p + integer, intent(in),dimension(2) :: pst + integer :: k + double precision, dimension(3) :: s + double precision, dimension(3) :: t + + ! Morse transform + do k=1,3 + t(k) = morse_transform(x(k), p, pst) + end do + s(1) = symmetrize(t(1), t(2), t(3), 'a') + s(2) = symmetrize(t(1), t(2), t(3), 'x') + s(3) = symmetrize(t(1), t(2), t(3), 'y') + end function morse_and_symmetrize + +! subroutine get_invariants(s, inv_out) +! use dim_parameter, only: qn +! use select_monom_mod, only: v_e_monom, v_ee_monom +! double precision, intent(in) :: s(qn) +! double precision, intent(out) :: inv_out(:) +! ! double precision, parameter :: ck = 1.00d0, dk = 1.00d0/ck ! scaling for higher order invariants +! double precision inv(24) +! integer, parameter :: inv_order(12) = & ! the order in which the invariants are selected +! & [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] +! double precision Rch, umb, xR, yR, xAng, yAng +!! for readability +! Rch = s(1) +! xR = s(2) +! yR = s(3) +! xAng = s(4) +! yAng = s(5) +! umb = s(6)**2 +!! invarianten +!! a moden +! inv(1) = Rch +! inv(2) = umb +!! invariante e pairs +! inv(3) = v_e_monom(xR, yR, 1) +! inv(4) = v_e_monom(xAng, yAng, 1) +!! third order e pairs +! inv(5) = v_e_monom(xR, yR, 2) +! inv(6) = v_e_monom(xAng, yAng, 2) +! ! invariant ee coupling +! inv(7) = v_ee_monom(xR, yR, xAng, yAng, 1) +! ! mode combinations +! inv(8) = Rch*umb +! +! inv(9) = Rch*v_e_monom(xR, yR, 1) +! inv(10) = umb*v_e_monom(xR, yR, 1) +! +! inv(11) = Rch*v_e_monom(xAng, yAng, 1) +! inv(12) = umb*v_e_monom(xAng, yAng, 1) +! +!! damp coordinates because of second order and higher invariants +! inv(3) = sign(sqrt(abs(inv(3))), inv(3)) +! inv(4) = sign(sqrt(abs(inv(4))), inv(4)) +! inv(5) = sign((abs(inv(5))**(1./3.)), inv(5)) +! inv(6) = sign((abs(inv(6))**(1./3.)), inv(6)) +! inv(7) = sign((abs(inv(7))**(1./3.)), inv(7)) +! inv(8) = sign(sqrt(abs(inv(8))), inv(8)) +! inv(9) = sign((abs(inv(9))**(1./3.)), inv(9)) +! inv(10) = sign((abs(inv(10))**(1./3.)), inv(10)) +! inv(11) = sign((abs(inv(11))**(1./3.)), inv(11)) +! inv(12) = sign((abs(inv(12))**(1./3.)), inv(12)) +! +! inv_out(:) = inv(inv_order(1:size(inv_out, 1))) +! +! end subroutine get_invariants + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % real part of spherical harmonics +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! Ylm shifted to 0 for theta=0 +! double precision function ylm(theta, phi, l, m) +! implicit none +! double precision theta, phi +! integer l, m +! ylm = plm2(dcos(theta), l, m)*cos(m*phi) - plm2(1.00d0, l, m) +! end function ylm +!!---------------------------------------------------------- +! double precision function plm2(x, l, n) +! implicit none +! double precision x +! integer l, m, n +! +! double precision pmm, p_mp1m, pllm +! integer ll +! +!! negative m und bereich von x abfangen +! if ((l .lt. 0)& +! &.or. (abs(n) .gt. abs(l))& +! &.or. (abs(x) .gt. 1.)) then +! write (6, '(''bad arguments in legendre'')') +! stop +! end if +! +!! fix sign of m to compute the positiv m +! m = abs(n) +! +! pmm = (-1)**m*dsqrt(fac(2*m))*1./((2**m)*fac(m))& !compute P(m,m) not P(l,l) +! &*(dsqrt(1.-x**2))**m +! +! if (l .eq. m) then +! plm2 = pmm !P(l,m)=P(m,m) +! else +! p_mp1m = x*dsqrt(dble(2*m + 1))*pmm !compute P(m+1,m) +! if (l .eq. m + 1) then +! plm2 = p_mp1m !P(l,m)=P(m+1,m) +! else +! do ll = m + 2, l +! pllm = x*(2*ll - 1)/dsqrt(dble(ll**2 - m**2))*p_mp1m& ! compute P(m+2,m) up to P(l,m) recursively +! &- dsqrt(dble((ll - 1)**2 - m**2))& +! &/dsqrt(dble(l**2 - m**2))*pmm +!! schreibe m+2 und m+1 jeweils fuer die naechste iteration +! pmm = p_mp1m !P(m,m) = P(m+1,m) +! p_mp1m = pllm !P(m+1,m) = P(m+2,m) +! end do +! plm2 = pllm !P(l,m)=P(m+k,m), k element N +! end if +! end if +! +!! sets the phase of -m term right (ignored to gurantee Ylm=(Yl-m)* for JT terms +!! if(n.lt.0) then +!! plm2 = (-1)**m * plm2 !* fac(l-m)/fac(l+m) +!! endif +! +! end function +!---------------------------------------------------------------------------------------------------- + double precision function fac(i) + integer i + select case (i) + case (0) + fac = 1.00d0 + case (1) + fac = 1.00d0 + case (2) + fac = 2.00d0 + case (3) + fac = 6.00d0 + case (4) + fac = 24.00d0 + case (5) + fac = 120.00d0 + case (6) + fac = 720.00d0 + case (7) + fac = 5040.00d0 + case (8) + fac = 40320.00d0 + case (9) + fac = 362880.00d0 + case (10) + fac = 3628800.00d0 + case (11) + fac = 39916800.00d0 + case (12) + fac = 479001600.00d0 + case default + write (*, *) 'ERROR: no case for given faculty, Max is 12!' + stop + end select + end function fac + + ! Does the simplest morse transform possible + ! one skaling factor + shift + function morse_transform(x, p, pst) result(t) + double precision, intent(in) :: x + double precision, intent(in) :: p(11) + integer, intent(in) :: pst(2) + double precision :: t + if (pst(2) == 11) then + t = 1.00d0 - exp(-abs(p(2))*(x - p(1))) + else + error stop 'in morse_transform key required or wrong number of parameters' + end if + end function morse_transform + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % FUNCTION F(...) ! MAIK DEPRICATING OVER THE TOP MORSE FUNCTION FOR MYSELF +! % +! % Returns exponent of tunable Morse coordinate +! % exponent is polynomial * gaussian (skewed) +! % ilabel = 1 or 2 selects the parameters a and sfac to be used +! % +! % Background: better representation of the prefector in the +! % exponend of the morse function. +! % Formular: f(r) = lest no3 paper +! % +! % Variables: +! % x: distance of atoms (double) +! % p: parameter vector (double[20]) +! % ii: 1 for CCl and 2 for CCH (int) +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + pure function f(x, p, ii) + integer , intent(in) :: ii !1 for CCL and 2 for CCH + double precision, intent(in) :: x !coordinate + double precision, intent(in) :: p(11) !parameter-vector + + integer i !running index + + double precision r !equilibrium distance + double precision gaus !gaus part of f + double precision poly !polynom part of f + double precision skew !tanh part of f + + double precision f !prefactor of exponent and returned value + + integer npoly(2) !order of polynom + +! Maximum polynom order + npoly(1) = 5 + npoly(2) = 5 + +! p(1): position of equilibrium +! p(2): constant of exponent +! p(3): constant for skewing the gaussian +! p(4): tuning for skewing the gaussian +! p(5): Gaussian exponent +! p(6): Shift of Gaussian maximum +! p(7)...: polynomial coefficients +! p(8+n)...: coefficients of Morse Power series + +! 1-exp{[p(2)+exp{-p(5)[x-p(6)]^2}[Taylor{p(7+n)}(x-p(6))]][x-p(1)]} + +! Tunable Morse function +! Power series in Tunable Morse coordinates of order m +! exponent is polynomial of order npoly * gaussian + switching function + +! set r r-r_e + r = x + r = r - p(1) + +! set up skewing function: + skew = 0.50d0*p(3)*(dtanh(dabs(p(4))*(r - p(6))) + 1.00d0) + +! set up gaussian function: + gaus = dexp(-dabs(p(5))*(r - p(6))**2) + +! set up power series: + poly = 0.00d0 + do i = 0, npoly(ii) - 1 + poly = poly + p(7 + i)*(r - p(6))**i + end do +! set up full exponent function: + f = dabs(p(2)) + skew + gaus*poly + + end function +!---------------------------------------------------------------------------------------------------- + pure function xproduct(a, b) result(axb) + double precision, intent(in) :: a(3), b(3) + double precision :: axb(3) !crossproduct a x b + axb(1) = a(2)*b(3) - a(3)*b(2) + axb(2) = a(3)*b(1) - a(1)*b(3) + axb(3) = a(1)*b(2) - a(2)*b(1) + end function xproduct + + pure function normalized(v) result(r) + double precision, intent(in) :: v(:) + double precision :: r(size(v)) + r = v/norm(v) + end function normalized + + pure function norm(v) result(n) + double precision, intent(in) :: v(:) + double precision n + n = dsqrt(sum(v(:)**2)) + end function norm + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % FUNCTION Project_Point_Into_Plane(x,n,r0) result(p) +! % return the to n orthogonal part of a vector x-r0 +! % p: projected point in plane +! % x: point being projected +! % n: normalvector of plane +! % r0: Point in plane +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + pure function project_point_into_plane(x, n, r0) result(p) + double precision, intent(in) :: x(:), n(:), r0(:) + double precision :: p(size(x)), xs(size(x)) + xs = x - r0 + p = xs - plane_to_point(x, n, r0) + end function project_point_into_plane + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % Function Plane_To_Point(x,n,r0) result(p) +! % p: part of n in x +! % x: point being projected +! % n: normalvector of plane +! % r0: Point in plane +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + pure function plane_to_point(x, n, r0) result(p) + double precision, intent(in) :: x(:), n(:), r0(:) + double precision p(size(x)), xs(size(x)), nn(size(n)) + nn = normalized(n) + xs = x - r0 + p = dot_product(nn, xs)*nn + end function plane_to_point + + subroutine check_coordinates(q) +! check for faulty kartesain coordinates + double precision, intent(in) :: q(:) + integer :: i + if (all(abs(q) <= epsilon(0.00d0))) then + stop 'Error (ctrans): all kartesian coordinates are<=1d-8' + end if + do i = 1, 9, 3 + if (all(abs(q(i:i + 2)) <= epsilon(0.00d0))) then + write (*, *) q + stop 'Error(ctrans):kartesian coordinates zero for one atom' + end if + end do + end subroutine + + pure function rotor_a_to_z(a, z) result(r) + double precision, intent(in) :: a(3), z(3) + double precision :: r(3, 3) + double precision :: alpha + double precision :: s1(3), s(3, 3), rotor(3, 3) + s1 = xproduct(normalized(a), normalized(z)) + alpha = asin(norm(s1)) + s(:, 1) = normalized(s1) + s(:, 2) = normalized(z) + s(:, 3) = xproduct(s1, z) + rotor = init_rotor(alpha, 0.00d0, 0.00d0) + r = matmul(s, matmul(rotor, transpose(s))) + end function + +! function returning Rz(gamma) * Ry(beta) * Rx(alpha) for basis order xyz + pure function init_rotor(alpha, beta, gamma) result(rotor) + double precision, intent(in) :: alpha, beta, gamma + double precision :: rotor(3, 3) + rotor = 0.00d0 + rotor(1, 1) = dcos(beta)*dcos(gamma) + rotor(1, 2) = dsin(alpha)*dsin(beta)*dcos(gamma)& + &- dcos(alpha)*dsin(gamma) + rotor(1, 3) = dcos(alpha)*dsin(beta)*dcos(gamma)& + &+ dsin(alpha)*dsin(gamma) + + rotor(2, 1) = dcos(beta)*dsin(gamma) + rotor(2, 2) = dsin(alpha)*dsin(beta)*dsin(gamma)& + &+ dcos(alpha)*dcos(gamma) + rotor(2, 3) = dcos(alpha)*dsin(beta)*dsin(gamma)& + &- dsin(alpha)*dcos(gamma) + + rotor(3, 1) = -dsin(beta) + rotor(3, 2) = dsin(alpha)*dcos(beta) + rotor(3, 3) = dcos(alpha)*dcos(beta) + end function init_rotor + + pure function create_plane(a, b, c) result(n) + double precision, intent(in) :: a(3), b(3), c(3) + double precision :: n(3) + double precision :: axb(3), bxc(3), cxa(3) + axb = xproduct(a, b) + bxc = xproduct(b, c) + cxa = xproduct(c, a) + n = normalized(axb + bxc + cxa) + end function create_plane + + function symmetrize(q1, q2, q3, sym) result(s) + double precision, intent(in) :: q1, q2, q3 + character, intent(in) :: sym + double precision :: s + select case (sym) + case ('a') + s = (q1 + q2 + q3)*sq3 + case ('x') + s = sq6*(2.00d0*q1 - q2 - q3) + case ('y') + s = sq2*(q2 - q3) + case default + write (*, *) 'ERROR: no rule for symmetrize with sym=', sym + stop + end select + end function symmetrize + + subroutine construct_HBend(ex, ey, ph1, ph2, ph3, x_axis, y_axis) + double precision, intent(in) :: ph1(3), ph2(3), ph3(3) + double precision, intent(in) :: x_axis(3), y_axis(3) + double precision, intent(out) :: ex, ey + double precision :: x1, y1, alpha1 + double precision :: x2, y2, alpha2 + double precision :: x3, y3, alpha3 +! get x and y components of projected points + x1 = dot_product(ph1, x_axis) + y1 = dot_product(ph1, y_axis) + x2 = dot_product(ph2, x_axis) + y2 = dot_product(ph2, y_axis) + x3 = dot_product(ph3, x_axis) + y3 = dot_product(ph3, y_axis) +! -> calculate H deformation angles + alpha3 = datan2(y2, x2) + alpha2 = -datan2(y3, x3) !-120*ang2rad +! write(*,*)' atan2' +! write(*,*) 'alpha2:' , alpha2/ang2rad +! write(*,*) 'alpha3:' , alpha3/ang2rad + if (alpha2 .lt. 0) alpha2 = alpha2 + pi2 + if (alpha3 .lt. 0) alpha3 = alpha3 + pi2 + alpha1 = (pi2 - alpha2 - alpha3) +! write(*,*)' fixed break line' +! write(*,*) 'alpha1:' , alpha1/ang2rad +! write(*,*) 'alpha2:' , alpha2/ang2rad +! write(*,*) 'alpha3:' , alpha3/ang2rad + alpha1 = alpha1 !- 120.00d0*ang2rad + alpha2 = alpha2 !- 120.00d0*ang2rad + alpha3 = alpha3 !- 120.00d0*ang2rad +! write(*,*)' delta alpha' +! write(*,*) 'alpha1:' , alpha1/ang2rad +! write(*,*) 'alpha2:' , alpha2/ang2rad +! write(*,*) 'alpha3:' , alpha3/ang2rad +! write(*,*) + +! construct symmetric and antisymmetric H angles + ex = symmetrize(alpha1, alpha2, alpha3, 'x') + ey = symmetrize(alpha1, alpha2, alpha3, 'y') + end subroutine construct_HBend + + pure function construct_umbrella(nh1, nh2, nh3, n)& + &result(umb) + double precision, intent(in) :: nh1(3), nh2(3), nh3(3) + double precision, intent(in) :: n(3) + double precision :: umb + double precision :: theta(3) +! calculate projections for umberella angle + theta(1) = dacos(dot_product(n, nh1)) + theta(2) = dacos(dot_product(n, nh2)) + theta(3) = dacos(dot_product(n, nh3)) +! construct umberella angle + umb = sum(theta(1:3))/3.00d0 - 90.00d0*ang2rad + end function construct_umbrella + + pure subroutine construct_sphericals& + &(theta, phi, cf, xaxis, yaxis, zaxis) + double precision, intent(in) :: cf(3), xaxis(3), yaxis(3), zaxis(3) + double precision, intent(out) :: theta, phi + double precision :: x, y, z, v(3) + v = normalized(cf) + x = dot_product(v, normalized(xaxis)) + y = dot_product(v, normalized(yaxis)) + z = dot_product(v, normalized(zaxis)) + theta = dacos(z) + phi = -datan2(y, x) + end subroutine construct_sphericals + +! subroutine int2kart(internal, kart) +! double precision, intent(in) :: internal(6) +! double precision, intent(out) :: kart(9) +! double precision :: h1x, h1y, h1z +! double precision :: h2x, h2y, h2z +! double precision :: h3x, h3y, h3z +! double precision :: dch0, dch1, dch2, dch3 +! double precision :: a1, a2, a3, wci +! +! kart = 0.00d0 +! dch1 = dchequi + sq3*internal(1) + 2*sq6*internal(2) +! dch2 = dchequi + sq3*internal(1) - sq6*internal(2) + sq2*internal(3) +! dch3 = dchequi + sq3*internal(1) - sq6*internal(2) - sq2*internal(3) +! a1 = 2*sq6*internal(4) +! a2 = -sq6*internal(4) + sq2*internal(5) +! a3 = -sq6*internal(4) - sq2*internal(5) +! wci = internal(6) +! +! ! Berechnung kartesische Koordinaten +! ! ----------------------- +! h1x = dch1*cos(wci*ang2rad) +! h1y = 0.0 +! h1z = -dch1*sin(wci*ang2rad) +! +! h3x = dch2*cos((a2 + 120)*ang2rad)*cos(wci*ang2rad) +! h3y = dch2*sin((a2 + 120)*ang2rad)*cos(wci*ang2rad) +! h3z = -dch2*sin(wci*ang2rad) +! +! h2x = dch3*cos((-a3 - 120)*ang2rad)*cos(wci*ang2rad) +! h2y = dch3*sin((-a3 - 120)*ang2rad)*cos(wci*ang2rad) +! h2z = -dch3*sin(wci*ang2rad) +! +! kart(1) = h1x +! kart(2) = h1y +! kart(3) = h1z +! kart(4) = h2x +! kart(5) = h2y +! kart(6) = h2z +! kart(7) = h3x +! kart(8) = h3y +! kart(9) = h3z +! end subroutine int2kart + +end module ctrans_mod +!**** Define coordinate transformation applied to the input before fit. +!*** +!*** +!****Conventions: +!*** +!*** ctrans: subroutine transforming a single point in coordinate space + + subroutine trans_in(pat_in,ntot) + use ctrans_mod, only: ctrans + implicit none + + include 'nnparams.incl' + !include 'nndbg.incl' + !include 'nncommon.incl' + + double precision pat_in(maxnin,maxpats) + integer ntot + + integer j + + do j=1,ntot + call ctrans(pat_in(:,j)) + ! FIRST ELEMENT OF PAT-IN ARE USED BY NEURON NETWORK + write(62,'(6f16.8)') pat_in(4:9,j) + enddo + end diff --git a/src/model/red_invariants.f90 b/src/model/red_invariants.f90 new file mode 100644 index 0000000..278051f --- /dev/null +++ b/src/model/red_invariants.f90 @@ -0,0 +1,72 @@ +Module invariants_mod + implicit none + contains + !---------------------------------------------------- + subroutine invariants(a,xs,ys,xb,yb,b,invar) + implicit none + double precision, intent(in) :: a, xs, ys, xb, yb, b + double precision, intent(out) :: invar(24) + complex(8) :: q1, q2 + LOGICAL,PARAMETER:: debg =.FALSE. + integer :: i + ! express the coordinate in complex + + q1 = dcmplx(xs, ys) + q2 = dcmplx(xb, yb) + + ! compute the invariants + invar(24) = a + invar(23) =b**2 + + ! INVARIANTS OF KIND II + !------------------------ + + invar(1) = dreal( q1 * conjg(q1) ) ! r11 + invar(2) = dreal( q1 * conjg(q2) ) ! r12 + invar(3) = dreal( q2 * conjg(q2) ) ! r22 + invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2 + + + !INVATIANTS OF KIND III + !------------------------ + + invar(5) = dreal( q1 * q1 * q1 ) ! r111 + invar(6) = dreal( q1 * q1 * q2 ) ! r112 + invar(7) = dreal( q1 * q2 * q2 ) ! r122 + invar(8) = dreal( q2 * q2 * q2 ) ! r222 + invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2 + invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2 + invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2 + invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222 + + ! INVARIANTS OF KIND IV + !------------------------- + + invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 )) + invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 )) + invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 )) + invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 )) + + ! INVARIANTS OF KIND V + !---------------------- + + invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 )) + invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 )) + invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 )) + invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 )) + invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 )) + invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 )) + + if (debg) then + write(*,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) + write(*,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) + write(*,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) + write(*,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) + + write(*,*)"THE INPUT COORDINATE IN COMPLEX REPRES" + write(*,*)"---------------------------------------" + write(*,*)"xs =",dreal(q1), "ys=",dimag(q1) + endif + + end subroutine invariants +end module invariants_mod