jnshuti
/
Gen_ANN_nh3_base
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compare commits

base: jnshuti:main
Branches Tags
jnshuti:main
jnshuti:Planar_nh3
jnshuti:Pyramidal_nh3
..
compare: jnshuti:Pyramidal_nh3
Branches Tags
jnshuti:Planar_nh3
jnshuti:Pyramidal_nh3
jnshuti:main

9 Commits
main ... Pyramidal_

Author SHA1 Message Date
jean paul nshuti 67961496d5 rename files 2025-10-10 11:28:25 +02:00
jean paul nshuti 30790a924d modify the filenames to differentiate planar and pyramidal geometry 2025-10-10 09:44:10 +02:00
jean paul nshuti 977dcab593 add the genetic parameter of pyramidal 2025-10-09 18:22:59 +02:00
jean paul nshuti b23f225ec8 change value of Wter 2025-10-09 17:40:33 +02:00
jean paul nshuti e747afbd40 modify the invariants for only few of them 2025-10-09 17:28:39 +02:00
jean paul nshuti 79ad45bb9d Merge branch 'main' into Planar_nh3 
I want the branch to have the same src code
 2025-10-09 14:30:49 +02:00
jean paul nshuti 7d03faaef1 modify the makefile for compilinng the planar 2025-10-09 13:03:42 +02:00
jean paul nshuti 6f1323639a fix the model and nnadia 2025-10-09 13:01:49 +02:00
jean paul nshuti ccf300a77e Add the model for planar geometry of NH3+ 2025-10-09 10:05:57 +02:00
10 changed files with 2413 additions and 8 deletions
Show Stats Expand all files Collapse all files

177
Makefile Normal file
View File

@ -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/
mdl = $(src)model/
# DEFINE THE PATH
VPATH = $(src)
VPATH += $(ldir)
VPATH += $(pdir)
VPATH += $(mdl)
# extra dir
extra_dirs = logs nnfits scans
gincl = -I$(src) -I$(pdir) -I$(ldir) -I$(mdl)
# 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_nh3_py.o nncoords_nh3_py.o genetic_param_nh3_pyr.o model_nh3_py.o nnadia_nh3.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

16
src/error.f
View File

@ -62,7 +62,7 @@
!-------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------
subroutine nnweight(wterr) subroutine nnweight(wterr,pat_out)
implicit none implicit none
! Evaluate system specific weighting for 1 pattern. ! Evaluate system specific weighting for 1 pattern.
@ -72,11 +72,11 @@
!include 'JTmod.incl' !include 'JTmod.incl'
double precision wterr(maxpout)!,pat_out(maxpout) double precision wterr(maxpout),pat_out(maxpout)
!double precision eref(nstat) !double precision eref(nstat)
!double precision wten(3) !double precision wten(3)
integer i
!integer j !integer j
!wten=1 !wten=1
@ -93,7 +93,15 @@
! weighting of energies ! weighting of energies
! wterr(j)=wdamp(pat_out(j)-eref(j))*wterr(j)*wten(j) ! wterr(j)=wdamp(pat_out(j)-eref(j))*wterr(j)*wten(j)
!enddo !enddo
wterr = 1.0d0 do i=1,size(wterr)
if (wterr(i) .eq. 0.0d0) then
wterr(i) =0.0d0
else
wterr(i) = 1.0d0
endif
enddo
pat_out=pat_out
!wterr = 1.0d0
contains contains
double precision function wdamp(dE) double precision function wdamp(dE)

7
src/mkNN.f
View File

@ -408,11 +408,10 @@
write(nnunit,newline) write(nnunit,newline)
do k=1,npat do k=1,npat
!call nnweight(wterr(1,k),pat_out(1,k)) ! JP no need for the weighting schemme used for energy
! right now wter is set to 1.0 in nnweights call nnweight(wterr(1,k),pat_out(1,k))
call nnweight(wterr(1,k))
enddo enddo
pat_out = pat_out ! to avoid unsed complain durinng compilation !pat_out = pat_out ! to avoid unsed complain durinng compilation
total=0 total=0
do k=1,sets do k=1,sets
write(nnunit,'(A10,I6.5)') '# Scan Nr.',k write(nnunit,'(A10,I6.5)') '# Scan Nr.',k

877
src/model/genetic_param_nh3_pyr.f90 Normal file
View File

@ -0,0 +1,877 @@
Module dip_param
IMPLICIT NONE
Integer,parameter :: np=66
Double precision :: p(66)
integer :: pst(2,400)
contains
SUBROUTINE init_dip_planar_data()
implicit none
p( 1)= 0.515449059D+00
p( 2)= 0.000000000D+00
p( 3)= 0.978476807D+00
p( 4)= 0.000000000D+00
p( 5)= -0.225153302D+01
p( 6)= 0.000000000D+00
p( 7)= -0.196655329D+00
p( 8)= 0.000000000D+00
p( 9)= 0.000000000D+00
p( 10)= 0.000000000D+00
p( 11)= 0.000000000D+00
p( 12)= 0.000000000D+00
p( 13)= -0.308320582D+01
p( 14)= 0.000000000D+00
p( 15)= 0.000000000D+00
p( 16)= 0.000000000D+00
p( 17)= 0.000000000D+00
p( 18)= 0.000000000D+00
p( 19)= 0.000000000D+00
p( 20)= 0.000000000D+00
p( 21)= 0.000000000D+00
p( 22)= 0.155576080D+00
p( 23)= 0.399570072D+00
p( 24)= 0.000000000D+00
p( 25)= -0.168532613D+00
p( 26)= 0.000000000D+00
p( 27)= 0.000000000D+00
p( 28)= 0.000000000D+00
p( 29)= 0.000000000D+00
p( 30)= 0.000000000D+00
p( 31)= 0.000000000D+00
p( 32)= 0.000000000D+00
p( 33)= 0.000000000D+00
p( 34)= 0.000000000D+00
p( 35)= 0.000000000D+00
p( 36)= 0.000000000D+00
p( 37)= 0.000000000D+00
p( 38)= 0.000000000D+00
p( 39)= 0.000000000D+00
p( 40)= 0.400000000D-01
p( 41)= -0.348763946D+00
p( 42)= 0.000000000D+00
p( 43)= -0.352820810D+00
p( 44)= 0.000000000D+00
p( 45)= 0.000000000D+00
p( 46)= -0.567363489D+00
p( 47)= -0.244801150D+00
p( 48)= 0.000000000D+00
p( 49)= 0.000000000D+00
p( 50)= 0.000000000D+00
p( 51)= -0.128708810D-01
p( 52)= 0.000000000D+00
p( 53)= 0.000000000D+00
p( 54)= 0.000000000D+00
p( 55)= -0.111374303D+01
p( 56)= 0.000000000D+00
p( 57)= -0.194029120D+01
p( 58)= 0.000000000D+00
p( 59)= 0.000000000D+00
p( 60)= 0.270000000D-01
p( 61)= -0.231405897D+00
p( 62)= 0.000000000D+00
p( 63)= -0.278734547D+00
p( 64)= 0.000000000D+00
p( 65)= 0.000000000D+00
p( 66)= -0.929361369D+00
pst(1, 1)= 1
pst(2, 1)= 2
pst(1, 2)= 3
pst(2, 2)= 2
pst(1, 3)= 5
pst(2, 3)= 2
pst(1, 4)= 7
pst(2, 4)= 3
pst(1, 5)= 10
pst(2, 5)= 3
pst(1, 6)= 13
pst(2, 6)= 3
pst(1, 7)= 16
pst(2, 7)= 2
pst(1, 8)= 18
pst(2, 8)= 2
pst(1, 9)= 20
pst(2, 9)= 2
pst(1, 10)= 22
pst(2, 10)= 1
pst(1, 11)= 23
pst(2, 11)= 2
pst(1, 12)= 25
pst(2, 12)= 5
pst(1, 13)= 30
pst(2, 13)= 10
pst(1, 14)= 40
pst(2, 14)= 1
pst(1, 15)= 41
pst(2, 15)= 2
pst(1, 16)= 43
pst(2, 16)= 4
pst(1, 17)= 47
pst(2, 17)= 8
pst(1, 18)= 55
pst(2, 18)= 2
pst(1, 19)= 57
pst(2, 19)= 3
pst(1, 20)= 60
pst(2, 20)= 1
pst(1, 21)= 61
pst(2, 21)= 2
pst(1, 22)= 63
pst(2, 22)= 4
pst(1, 23)= 0
pst(2, 23)= 0
pst(1, 24)= 0
pst(2, 24)= 0
pst(1, 25)= 0
pst(2, 25)= 0
pst(1, 26)= 0
pst(2, 26)= 0
pst(1, 27)= 0
pst(2, 27)= 0
pst(1, 28)= 0
pst(2, 28)= 0
pst(1, 29)= 0
pst(2, 29)= 0
pst(1, 30)= 0
pst(2, 30)= 0
pst(1, 31)= 0
pst(2, 31)= 0
pst(1, 32)= 0
pst(2, 32)= 0
pst(1, 33)= 0
pst(2, 33)= 0
pst(1, 34)= 0
pst(2, 34)= 0
pst(1, 35)= 0
pst(2, 35)= 0
pst(1, 36)= 0
pst(2, 36)= 0
pst(1, 37)= 0
pst(2, 37)= 0
pst(1, 38)= 0
pst(2, 38)= 0
pst(1, 39)= 0
pst(2, 39)= 0
pst(1, 40)= 0
pst(2, 40)= 0
pst(1, 41)= 0
pst(2, 41)= 0
pst(1, 42)= 0
pst(2, 42)= 0
pst(1, 43)= 0
pst(2, 43)= 0
pst(1, 44)= 0
pst(2, 44)= 0
pst(1, 45)= 0
pst(2, 45)= 0
pst(1, 46)= 0
pst(2, 46)= 0
pst(1, 47)= 0
pst(2, 47)= 0
pst(1, 48)= 0
pst(2, 48)= 0
pst(1, 49)= 0
pst(2, 49)= 0
pst(1, 50)= 0
pst(2, 50)= 0
pst(1, 51)= 0
pst(2, 51)= 0
pst(1, 52)= 0
pst(2, 52)= 0
pst(1, 53)= 0
pst(2, 53)= 0
pst(1, 54)= 0
pst(2, 54)= 0
pst(1, 55)= 0
pst(2, 55)= 0
pst(1, 56)= 0
pst(2, 56)= 0
pst(1, 57)= 0
pst(2, 57)= 0
pst(1, 58)= 0
pst(2, 58)= 0
pst(1, 59)= 0
pst(2, 59)= 0
pst(1, 60)= 0
pst(2, 60)= 0
pst(1, 61)= 0
pst(2, 61)= 0
pst(1, 62)= 0
pst(2, 62)= 0
pst(1, 63)= 0
pst(2, 63)= 0
pst(1, 64)= 0
pst(2, 64)= 0
pst(1, 65)= 0
pst(2, 65)= 0
pst(1, 66)= 0
pst(2, 66)= 0
pst(1, 67)= 0
pst(2, 67)= 0
pst(1, 68)= 0
pst(2, 68)= 0
pst(1, 69)= 0
pst(2, 69)= 0
pst(1, 70)= 0
pst(2, 70)= 0
pst(1, 71)= 0
pst(2, 71)= 0
pst(1, 72)= 0
pst(2, 72)= 0
pst(1, 73)= 0
pst(2, 73)= 0
pst(1, 74)= 0
pst(2, 74)= 0
pst(1, 75)= 0
pst(2, 75)= 0
pst(1, 76)= 0
pst(2, 76)= 0
pst(1, 77)= 0
pst(2, 77)= 0
pst(1, 78)= 0
pst(2, 78)= 0
pst(1, 79)= 0
pst(2, 79)= 0
pst(1, 80)= 0
pst(2, 80)= 0
pst(1, 81)= 0
pst(2, 81)= 0
pst(1, 82)= 0
pst(2, 82)= 0
pst(1, 83)= 0
pst(2, 83)= 0
pst(1, 84)= 0
pst(2, 84)= 0
pst(1, 85)= 0
pst(2, 85)= 0
pst(1, 86)= 0
pst(2, 86)= 0
pst(1, 87)= 0
pst(2, 87)= 0
pst(1, 88)= 0
pst(2, 88)= 0
pst(1, 89)= 0
pst(2, 89)= 0
pst(1, 90)= 0
pst(2, 90)= 0
pst(1, 91)= 0
pst(2, 91)= 0
pst(1, 92)= 0
pst(2, 92)= 0
pst(1, 93)= 0
pst(2, 93)= 0
pst(1, 94)= 0
pst(2, 94)= 0
pst(1, 95)= 0
pst(2, 95)= 0
pst(1, 96)= 0
pst(2, 96)= 0
pst(1, 97)= 0
pst(2, 97)= 0
pst(1, 98)= 0
pst(2, 98)= 0
pst(1, 99)= 0
pst(2, 99)= 0
pst(1,100)= 0
pst(2,100)= 0
pst(1,101)= 0
pst(2,101)= 0
pst(1,102)= 0
pst(2,102)= 0
pst(1,103)= 0
pst(2,103)= 0
pst(1,104)= 0
pst(2,104)= 0
pst(1,105)= 0
pst(2,105)= 0
pst(1,106)= 0
pst(2,106)= 0
pst(1,107)= 0
pst(2,107)= 0
pst(1,108)= 0
pst(2,108)= 0
pst(1,109)= 0
pst(2,109)= 0
pst(1,110)= 0
pst(2,110)= 0
pst(1,111)= 0
pst(2,111)= 0
pst(1,112)= 0
pst(2,112)= 0
pst(1,113)= 0
pst(2,113)= 0
pst(1,114)= 0
pst(2,114)= 0
pst(1,115)= 0
pst(2,115)= 0
pst(1,116)= 0
pst(2,116)= 0
pst(1,117)= 0
pst(2,117)= 0
pst(1,118)= 0
pst(2,118)= 0
pst(1,119)= 0
pst(2,119)= 0
pst(1,120)= 0
pst(2,120)= 0
pst(1,121)= 0
pst(2,121)= 0
pst(1,122)= 0
pst(2,122)= 0
pst(1,123)= 0
pst(2,123)= 0
pst(1,124)= 0
pst(2,124)= 0
pst(1,125)= 0
pst(2,125)= 0
pst(1,126)= 0
pst(2,126)= 0
pst(1,127)= 0
pst(2,127)= 0
pst(1,128)= 0
pst(2,128)= 0
pst(1,129)= 0
pst(2,129)= 0
pst(1,130)= 0
pst(2,130)= 0
pst(1,131)= 0
pst(2,131)= 0
pst(1,132)= 0
pst(2,132)= 0
pst(1,133)= 0
pst(2,133)= 0
pst(1,134)= 0
pst(2,134)= 0
pst(1,135)= 0
pst(2,135)= 0
pst(1,136)= 0
pst(2,136)= 0
pst(1,137)= 0
pst(2,137)= 0
pst(1,138)= 0
pst(2,138)= 0
pst(1,139)= 0
pst(2,139)= 0
pst(1,140)= 0
pst(2,140)= 0
pst(1,141)= 0
pst(2,141)= 0
pst(1,142)= 0
pst(2,142)= 0
pst(1,143)= 0
pst(2,143)= 0
pst(1,144)= 0
pst(2,144)= 0
pst(1,145)= 0
pst(2,145)= 0
pst(1,146)= 0
pst(2,146)= 0
pst(1,147)= 0
pst(2,147)= 0
pst(1,148)= 0
pst(2,148)= 0
pst(1,149)= 0
pst(2,149)= 0
pst(1,150)= 0
pst(2,150)= 0
pst(1,151)= 0
pst(2,151)= 0
pst(1,152)= 0
pst(2,152)= 0
pst(1,153)= 0
pst(2,153)= 0
pst(1,154)= 0
pst(2,154)= 0
pst(1,155)= 0
pst(2,155)= 0
pst(1,156)= 0
pst(2,156)= 0
pst(1,157)= 0
pst(2,157)= 0
pst(1,158)= 0
pst(2,158)= 0
pst(1,159)= 0
pst(2,159)= 0
pst(1,160)= 0
pst(2,160)= 0
pst(1,161)= 0
pst(2,161)= 0
pst(1,162)= 0
pst(2,162)= 0
pst(1,163)= 0
pst(2,163)= 0
pst(1,164)= 0
pst(2,164)= 0
pst(1,165)= 0
pst(2,165)= 0
pst(1,166)= 0
pst(2,166)= 0
pst(1,167)= 0
pst(2,167)= 0
pst(1,168)= 0
pst(2,168)= 0
pst(1,169)= 0
pst(2,169)= 0
pst(1,170)= 0
pst(2,170)= 0
pst(1,171)= 0
pst(2,171)= 0
pst(1,172)= 0
pst(2,172)= 0
pst(1,173)= 0
pst(2,173)= 0
pst(1,174)= 0
pst(2,174)= 0
pst(1,175)= 0
pst(2,175)= 0
pst(1,176)= 0
pst(2,176)= 0
pst(1,177)= 0
pst(2,177)= 0
pst(1,178)= 0
pst(2,178)= 0
pst(1,179)= 0
pst(2,179)= 0
pst(1,180)= 0
pst(2,180)= 0
pst(1,181)= 0
pst(2,181)= 0
pst(1,182)= 0
pst(2,182)= 0
pst(1,183)= 0
pst(2,183)= 0
pst(1,184)= 0
pst(2,184)= 0
pst(1,185)= 0
pst(2,185)= 0
pst(1,186)= 0
pst(2,186)= 0
pst(1,187)= 0
pst(2,187)= 0
pst(1,188)= 0
pst(2,188)= 0
pst(1,189)= 0
pst(2,189)= 0
pst(1,190)= 0
pst(2,190)= 0
pst(1,191)= 0
pst(2,191)= 0
pst(1,192)= 0
pst(2,192)= 0
pst(1,193)= 0
pst(2,193)= 0
pst(1,194)= 0
pst(2,194)= 0
pst(1,195)= 0
pst(2,195)= 0
pst(1,196)= 0
pst(2,196)= 0
pst(1,197)= 0
pst(2,197)= 0
pst(1,198)= 0
pst(2,198)= 0
pst(1,199)= 0
pst(2,199)= 0
pst(1,200)= 0
pst(2,200)= 0
pst(1,201)= 0
pst(2,201)= 0
pst(1,202)= 0
pst(2,202)= 0
pst(1,203)= 0
pst(2,203)= 0
pst(1,204)= 0
pst(2,204)= 0
pst(1,205)= 0
pst(2,205)= 0
pst(1,206)= 0
pst(2,206)= 0
pst(1,207)= 0
pst(2,207)= 0
pst(1,208)= 0
pst(2,208)= 0
pst(1,209)= 0
pst(2,209)= 0
pst(1,210)= 0
pst(2,210)= 0
pst(1,211)= 0
pst(2,211)= 0
pst(1,212)= 0
pst(2,212)= 0
pst(1,213)= 0
pst(2,213)= 0
pst(1,214)= 0
pst(2,214)= 0
pst(1,215)= 0
pst(2,215)= 0
pst(1,216)= 0
pst(2,216)= 0
pst(1,217)= 0
pst(2,217)= 0
pst(1,218)= 0
pst(2,218)= 0
pst(1,219)= 0
pst(2,219)= 0
pst(1,220)= 0
pst(2,220)= 0
pst(1,221)= 0
pst(2,221)= 0
pst(1,222)= 0
pst(2,222)= 0
pst(1,223)= 0
pst(2,223)= 0
pst(1,224)= 0
pst(2,224)= 0
pst(1,225)= 0
pst(2,225)= 0
pst(1,226)= 0
pst(2,226)= 0
pst(1,227)= 0
pst(2,227)= 0
pst(1,228)= 0
pst(2,228)= 0
pst(1,229)= 0
pst(2,229)= 0
pst(1,230)= 0
pst(2,230)= 0
pst(1,231)= 0
pst(2,231)= 0
pst(1,232)= 0
pst(2,232)= 0
pst(1,233)= 0
pst(2,233)= 0
pst(1,234)= 0
pst(2,234)= 0
pst(1,235)= 0
pst(2,235)= 0
pst(1,236)= 0
pst(2,236)= 0
pst(1,237)= 0
pst(2,237)= 0
pst(1,238)= 0
pst(2,238)= 0
pst(1,239)= 0
pst(2,239)= 0
pst(1,240)= 0
pst(2,240)= 0
pst(1,241)= 0
pst(2,241)= 0
pst(1,242)= 0
pst(2,242)= 0
pst(1,243)= 0
pst(2,243)= 0
pst(1,244)= 0
pst(2,244)= 0
pst(1,245)= 0
pst(2,245)= 0
pst(1,246)= 0
pst(2,246)= 0
pst(1,247)= 0
pst(2,247)= 0
pst(1,248)= 0
pst(2,248)= 0
pst(1,249)= 0
pst(2,249)= 0
pst(1,250)= 0
pst(2,250)= 0
pst(1,251)= 0
pst(2,251)= 0
pst(1,252)= 0
pst(2,252)= 0
pst(1,253)= 0
pst(2,253)= 0
pst(1,254)= 0
pst(2,254)= 0
pst(1,255)= 0
pst(2,255)= 0
pst(1,256)= 0
pst(2,256)= 0
pst(1,257)= 0
pst(2,257)= 0
pst(1,258)= 0
pst(2,258)= 0
pst(1,259)= 0
pst(2,259)= 0
pst(1,260)= 0
pst(2,260)= 0
pst(1,261)= 0
pst(2,261)= 0
pst(1,262)= 0
pst(2,262)= 0
pst(1,263)= 0
pst(2,263)= 0
pst(1,264)= 0
pst(2,264)= 0
pst(1,265)= 0
pst(2,265)= 0
pst(1,266)= 0
pst(2,266)= 0
pst(1,267)= 0
pst(2,267)= 0
pst(1,268)= 0
pst(2,268)= 0
pst(1,269)= 0
pst(2,269)= 0
pst(1,270)= 0
pst(2,270)= 0
pst(1,271)= 0
pst(2,271)= 0
pst(1,272)= 0
pst(2,272)= 0
pst(1,273)= 0
pst(2,273)= 0
pst(1,274)= 0
pst(2,274)= 0
pst(1,275)= 0
pst(2,275)= 0
pst(1,276)= 0
pst(2,276)= 0
pst(1,277)= 0
pst(2,277)= 0
pst(1,278)= 0
pst(2,278)= 0
pst(1,279)= 0
pst(2,279)= 0
pst(1,280)= 0
pst(2,280)= 0
pst(1,281)= 0
pst(2,281)= 0
pst(1,282)= 0
pst(2,282)= 0
pst(1,283)= 0
pst(2,283)= 0
pst(1,284)= 0
pst(2,284)= 0
pst(1,285)= 0
pst(2,285)= 0
pst(1,286)= 0
pst(2,286)= 0
pst(1,287)= 0
pst(2,287)= 0
pst(1,288)= 0
pst(2,288)= 0
pst(1,289)= 0
pst(2,289)= 0
pst(1,290)= 0
pst(2,290)= 0
pst(1,291)= 0
pst(2,291)= 0
pst(1,292)= 0
pst(2,292)= 0
pst(1,293)= 0
pst(2,293)= 0
pst(1,294)= 0
pst(2,294)= 0
pst(1,295)= 0
pst(2,295)= 0
pst(1,296)= 0
pst(2,296)= 0
pst(1,297)= 0
pst(2,297)= 0
pst(1,298)= 0
pst(2,298)= 0
pst(1,299)= 0
pst(2,299)= 0
pst(1,300)= 0
pst(2,300)= 0
pst(1,301)= 0
pst(2,301)= 0
pst(1,302)= 0
pst(2,302)= 0
pst(1,303)= 0
pst(2,303)= 0
pst(1,304)= 0
pst(2,304)= 0
pst(1,305)= 0
pst(2,305)= 0
pst(1,306)= 0
pst(2,306)= 0
pst(1,307)= 0
pst(2,307)= 0
pst(1,308)= 0
pst(2,308)= 0
pst(1,309)= 0
pst(2,309)= 0
pst(1,310)= 0
pst(2,310)= 0
pst(1,311)= 0
pst(2,311)= 0
pst(1,312)= 0
pst(2,312)= 0
pst(1,313)= 0
pst(2,313)= 0
pst(1,314)= 0
pst(2,314)= 0
pst(1,315)= 0
pst(2,315)= 0
pst(1,316)= 0
pst(2,316)= 0
pst(1,317)= 0
pst(2,317)= 0
pst(1,318)= 0
pst(2,318)= 0
pst(1,319)= 0
pst(2,319)= 0
pst(1,320)= 0
pst(2,320)= 0
pst(1,321)= 0
pst(2,321)= 0
pst(1,322)= 0
pst(2,322)= 0
pst(1,323)= 0
pst(2,323)= 0
pst(1,324)= 0
pst(2,324)= 0
pst(1,325)= 0
pst(2,325)= 0
pst(1,326)= 0
pst(2,326)= 0
pst(1,327)= 0
pst(2,327)= 0
pst(1,328)= 0
pst(2,328)= 0
pst(1,329)= 0
pst(2,329)= 0
pst(1,330)= 0
pst(2,330)= 0
pst(1,331)= 0
pst(2,331)= 0
pst(1,332)= 0
pst(2,332)= 0
pst(1,333)= 0
pst(2,333)= 0
pst(1,334)= 0
pst(2,334)= 0
pst(1,335)= 0
pst(2,335)= 0
pst(1,336)= 0
pst(2,336)= 0
pst(1,337)= 0
pst(2,337)= 0
pst(1,338)= 0
pst(2,338)= 0
pst(1,339)= 0
pst(2,339)= 0
pst(1,340)= 0
pst(2,340)= 0
pst(1,341)= 0
pst(2,341)= 0
pst(1,342)= 0
pst(2,342)= 0
pst(1,343)= 0
pst(2,343)= 0
pst(1,344)= 0
pst(2,344)= 0
pst(1,345)= 0
pst(2,345)= 0
pst(1,346)= 0
pst(2,346)= 0
pst(1,347)= 0
pst(2,347)= 0
pst(1,348)= 0
pst(2,348)= 0
pst(1,349)= 0
pst(2,349)= 0
pst(1,350)= 0
pst(2,350)= 0
pst(1,351)= 0
pst(2,351)= 0
pst(1,352)= 0
pst(2,352)= 0
pst(1,353)= 0
pst(2,353)= 0
pst(1,354)= 0
pst(2,354)= 0
pst(1,355)= 0
pst(2,355)= 0
pst(1,356)= 0
pst(2,356)= 0
pst(1,357)= 0
pst(2,357)= 0
pst(1,358)= 0
pst(2,358)= 0
pst(1,359)= 0
pst(2,359)= 0
pst(1,360)= 0
pst(2,360)= 0
pst(1,361)= 0
pst(2,361)= 0
pst(1,362)= 0
pst(2,362)= 0
pst(1,363)= 0
pst(2,363)= 0
pst(1,364)= 0
pst(2,364)= 0
pst(1,365)= 0
pst(2,365)= 0
pst(1,366)= 0
pst(2,366)= 0
pst(1,367)= 0
pst(2,367)= 0
pst(1,368)= 0
pst(2,368)= 0
pst(1,369)= 0
pst(2,369)= 0
pst(1,370)= 0
pst(2,370)= 0
pst(1,371)= 0
pst(2,371)= 0
pst(1,372)= 0
pst(2,372)= 0
pst(1,373)= 0
pst(2,373)= 0
pst(1,374)= 0
pst(2,374)= 0
pst(1,375)= 0
pst(2,375)= 0
pst(1,376)= 0
pst(2,376)= 0
pst(1,377)= 0
pst(2,377)= 0
pst(1,378)= 0
pst(2,378)= 0
pst(1,379)= 0
pst(2,379)= 0
pst(1,380)= 0
pst(2,380)= 0
pst(1,381)= 0
pst(2,381)= 0
pst(1,382)= 0
pst(2,382)= 0
pst(1,383)= 0
pst(2,383)= 0
pst(1,384)= 0
pst(2,384)= 0
pst(1,385)= 0
pst(2,385)= 0
pst(1,386)= 0
pst(2,386)= 0
pst(1,387)= 0
pst(2,387)= 0
pst(1,388)= 0
pst(2,388)= 0
pst(1,389)= 0
pst(2,389)= 0
pst(1,390)= 0
pst(2,390)= 0
pst(1,391)= 0
pst(2,391)= 0
pst(1,392)= 0
pst(2,392)= 0
pst(1,393)= 0
pst(2,393)= 0
pst(1,394)= 0
pst(2,394)= 0
pst(1,395)= 0
pst(2,395)= 0
pst(1,396)= 0
pst(2,396)= 0
pst(1,397)= 0
pst(2,397)= 0
pst(1,398)= 0
pst(2,398)= 0
pst(1,399)= 0
pst(2,399)= 0
pst(1,400)= 0
pst(2,400)= 0
End SUBROUTINE init_dip_planar_data
End Module dip_param

82
src/model/invariants_nh3_py.f90 Normal file
View File

@ -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(4)
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(1)
inv(2) = invar(5)
inv(3) = invar(9)
inv(4) = invar(23)
if (debg) then
write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4)
write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12)
write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16)
write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22)
write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES"
write(14,*)"---------------------------------------"
write(14,*)"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

405
src/model/model_nh3_py.f90 Normal file
View File

@ -0,0 +1,405 @@
module diabmodel
use iso_fortran_env, only: idp => int32, dp => real64
use dip_param
implicit none
include "nnparams.incl"
integer(idp),parameter:: ndiab=4
logical :: debug=.false.
contains
subroutine diab_x(e,q,nn_out)
real(dp),intent(in)::q(maxnin)
real(dp),intent(inout):: nn_out(maxnout)
real(dp),intent(out)::e(ndiab,ndiab)
integer(idp) id,i,j ,ii
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
real(dp),dimension(21):: shift,scal
xs=q(6)
ys=q(7)
xb=q(8)
yb=q(9)
a=q(5)
b=q(10)
call init_dip_planar_data()
! modify the parametr
!shift=1.0_dp
!scal=1.0d-2
scal = [81.50814,1.19850,245.66937,885.73515,875.91907,4.55261,4.14318,17.98337, &
47.60657,24.68377,124.20271,40.13721,153.81484,13.19697,258.11015,769.62299, &
147.59778,97.94602,413.29279,62.33428,1.88373]
shift = [-45.54845,-0.12528,139.98528,887.54027,602.90180,-1.61771,2.28343,-12.43082, &
-29.81680,9.29555,123.83516,-16.80044,108.07947,-2.19734,120.51791,496.11692, &
-83.20774,43.09530,373.00489,-15.32774,1.30491]
! V term of A2''
ii=1
do i =1,np
if (p(i) .ne. 0.0d0) then
p(i) =p(i)*(shift(ii) + scal(ii)*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
v3_vec( 1) = xs*(xs**2-3*ys**2)
v3_vec( 2) = xb*(xb**2-3*yb**2)
v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
v3_vec( 5) = ys*(3*xs**2-ys**2)
v3_vec( 6) = yb*(3*xb**2-yb**2)
v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb
e=0.0d0
id=1 !1
! V-term
! order 1
e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
id=id+1 !2
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
id=id+1 !3
e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
! order 2
id=id+1 !4
e(1,1)=e(1,1)+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 !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)
! 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
id=id+1 !8
e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
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
! JAHN TELLER COUPLING W AND Z
! 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))
! order 1
id=id+1 !11
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
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)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb
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)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb)
e(2,3)=e(2,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)
! order 3
id=id+1 !13
do i=1,4
j=i-1
e(2,2)=e(2,2)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
e(3,3)=e(3,3)-(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
e(2,3)=e(2,3)+(-p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i+4)
enddo
e(2,2)=e(2,2)+p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb
e(3,3)=e(3,3)-(p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb)
e(2,3)=e(2,3)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb
! PSEUDO JAHN TELLER
! A2 ground state coupled with E
! ###################################################
! ###################################################
! order 0
id=id+1 !14
e(1,2)=e(1,2)+b*p(pst(1,id))
! order 1
id=id+1 !15
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 !16
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) + p(pst(1,id)+3)*(xs**2+ys**2))
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))
! order 3
id =id+1 ! 17
do i=1,4
e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
e(1,3)=e(1,3)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
enddo
!! THE COUPLING OF A2 WITH A1
!####################################################
!####################################################
! order 1
id=id+1 !18
e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
id=id+1 !19
e(1,4)=e(1,4)+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))
!!! THE COUPLING OF A1 WITH E
!!####################################################
!####################################################
! order 0
id=id+1 !20
e(2,4)=e(2,4)+p(pst(1,id))
! order 1
id=id+1 !21
e(2,4)=e(2,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
e(3,4)=e(3,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
! order 2
id=id+1 !22
e(2,4)=e(2,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) +p(pst(1,id)+3)*(xs**2+ys**2)
e(3,4)=e(3,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)
! order 3
id=id+1 !23
do i=1,4
e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
e(3,4)=e(3,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
enddo
!! End of the model
e(2,1)=e(1,2)
e(3,1)=e(1,3)
e(3,2)=e(2,3)
e(4,1)=e(1,4)
e(4,2)=e(2,4)
e(4,3)=e(3,4)
end subroutine diab_x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! THE Y COMPONENT OF DIPOLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine diab_y(e,q,nn_out)
real(dp),intent(in)::q(maxnin)
real(dp),intent(inout):: nn_out(maxnout)
real(dp),intent(out)::e(ndiab,ndiab)
integer(idp) id,i,j, ii
real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
real(dp),dimension(21):: shift,scal
xs=q(6)
ys=q(7)
xb=q(8)
yb=q(9)
a=q(5)
b=q(10)
call init_dip_planar_data()
! modify the parametr
!shift=1.0_dp
!scal=1.0d-2
scal = [81.50814,1.19850,245.66937,885.73515,875.91907,4.55261,4.14318,17.98337, &
47.60657,24.68377,124.20271,40.13721,153.81484,13.19697,258.11015,769.62299, &
147.59778,97.94602,413.29279,62.33428,1.88373]
shift = [-45.54845,-0.12528,139.98528,887.54027,602.90180,-1.61771,2.28343,-12.43082, &
-29.81680,9.29555,123.83516,-16.80044,108.07947,-2.19734,120.51791,496.11692, &
-83.20774,43.09530,373.00489,-15.32774,1.30491]
! V term of A2''
ii=1
do i =1,np
if (p(i) .ne. 0) then
p(i) =p(i)*(shift(ii) + scal(ii)*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
v3_vec( 1) = xs*(xs**2-3*ys**2)
v3_vec( 2) = xb*(xb**2-3*yb**2)
v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
v3_vec( 5) = ys*(3*xs**2-ys**2)
v3_vec( 6) = yb*(3*xb**2-yb**2)
v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb
e=0.0d0
! V-term
id=1 !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
! order 2
id=id+1 !4
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)
! 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 !8
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 !9
e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb
! V- term + totally symmetric coord a
! JAHN TELLER COUPLING TERM
! order 0
id=id+1 !10
e(2,3)=e(2,3)+p(pst(1,id))
! order 1
id=id+1 !11
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
!id=id+1 !12
! order 2
id=id+1 !12
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)
e(2,3)=e(2,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)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb
! order 3
id=id+1 !13
do i=1,4
j=i-1
e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4)
e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4)
e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
enddo
e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb
e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb
e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb
! PSEUDO JAHN TELLER
! ORDER 0
! THE COUPLING OF A2 GROUND STATE WITH E
! ###################################################
! ###################################################
! order 0
id=id+1 !14
e(1,3)=e(1,3)-b*(p(pst(1,id)))
! order 1
id=id+1 !15
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 !16
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) - p(pst(1,id)+3)*(xs**2+ys**2))
! order 3
id = id+1 ! 17
do i=1,4
e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
e(1,3)=e(1,3)-b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
enddo
! THE COUPLING OF A2 WITH A1
!####################################################
!####################################################
! order 1
id=id+1 !17
e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
! order 2
id=id+1 !18
e(1,4)=e(1,4)-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))
! THE COUPLING OF A1 WITH E
!####################################################
!####################################################
! order 0
id=id+1 !19
e(3,4)=e(3,4)-p(pst(1,id))
! order 1
id=id+1 !20
e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb
e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
! order 2
id=id+1 !21
e(2,4)=e(2,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(3,4)=e(3,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) - p(pst(1,id)+3)*(xs**2+ys**2)
id =id+1 ! 23
! order 3
do i=1,4
e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
e(3,4)=e(3,4)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
enddo
! end of the model
e(2,1)=e(1,2)
e(3,1)=e(1,3)
e(3,2)=e(2,3)
e(4,1)=e(1,4)
e(4,2)=e(2,4)
e(4,3)=e(3,4)
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 = 2, size(mat, 1)
do m = 1, n - 1
mat(n, m) = mat(m, n)
end do
end do
end subroutine copy_2_lower_triangle
end module diabmodel

68
src/model/nnadia_nh3.f Normal file
View File

@ -0,0 +1,68 @@
subroutine nninit_mod()
implicit none
include 'nnparams.incl'
end subroutine
subroutine nnadia(coords,coeffs,adiaoutp)
USE diabmodel, only: diab_x, diab_y
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=4
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(dmat_x,coords,coeffs)
call diab_y(dmat_y,coords,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(1,2)
adiaoutp(3) = dmat_x(1,3)
adiaoutp(4) = dmat_x(1,4)
adiaoutp(5) = dmat_x(2,2)
adiaoutp(6) = dmat_x(2,3)
adiaoutp(7) = dmat_x(2,4)
adiaoutp(8) = dmat_x(3,3)
adiaoutp(9) = dmat_x(3,4)
adiaoutp(10) = dmat_x(4,4)
adiaoutp(11) = -1.0*dmat_y(1,1)
adiaoutp(12) = -1.0*dmat_y(1,2)
adiaoutp(13) = -1.0*dmat_y(1,3)
adiaoutp(14) = -1.0*dmat_y(1,4)
adiaoutp(15) = -1.0*dmat_y(2,2)
adiaoutp(16) = -1.0*dmat_y(2,3)
adiaoutp(17) = -1.0*dmat_y(2,4)
adiaoutp(18) = -1.0*dmat_y(3,3)
adiaoutp(19) = -1.0*dmat_y(3,4)
adiaoutp(20) = -1.0*dmat_y(4,4)
!write(*,*) dmat_x(1,1)
END SUBROUTINE

716
src/model/nncoords_nh3_py.f90 Normal file
View File

@ -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 = 1.0228710942d0 ! req NH3+
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(4)
! debugging
logical, parameter :: dbg = .false.
! initialize coordinate vectors
s = 0.0d0
t = 0.0d0
! write kartesian coords for readability
c_atom(1:3) = 0.0d0 ! N-atom at origin
do k = 1, 3
ch1(k) = q(k )
ch2(k) = q(k + 3)
ch3(k) = q(k + 6)
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(aR,exR,eyR,exAng,eyAng,umb,inv)
q(1:4)=inv(1:4)
q(5) = aR
q(6) = exR
q(7) = eyR
q(8) = exAng
q(9) = eyAng
q(10) = -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

1
src/model/nnparams.incl Symbolic link
View File

@ -0,0 +1 @@
../nnparams.incl

72
src/model/red_invariants.f90 Normal file
View File

@ -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