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

Compare commits

base: jnshuti:main
Branches Tags
jnshuti:main
jnshuti:model_muxy_nh3
jnshuti:model_lmat_no3
..
compare: jnshuti:model_lmat_no3
Branches Tags
jnshuti:model_muxy_nh3
jnshuti:main
jnshuti:model_lmat_no3

1 Commits
main ... model_lmat

Author SHA1 Message Date
jean paul nshuti de65a49694 first commit of the branch of lmatrix of no3 2025-10-08 15:23:56 +02:00
13 changed files with 2756 additions and 0 deletions
Show Stats Expand all files Collapse all files

173
Makefile Normal file
View File

@ -0,0 +1,173 @@
SHELL = /bin/bash
.SUFFIXES :
.SUFFIXES : .f .o .f90
src = ./src/
build = ./obj/
bin = ./bin/
######################################################################
version=localdw-1.0
######################################################################
#IFORT VERSION (DEFAULT)
FC = ifort
#MODERN IFORT VERSION (for compiling on laptops)
FFLAGS =-O2 -qopenmp -qmkl -heap-arrays -module $(build) -cpp -g -diag-disable=10448
#-openmp -complex-limited-range -xW -i-static -ip -ftz -no-prec-div -opt-prefetch -heap-arrays -align dcommons -mkl -mcmodel=medium
DBGFLAGS = -debug -check -check bounds #-warn uncalled -warn nousage -warn nounused -openmp -warn -warn notruncated_source
DBGFLAGS+= -pg
#GFORTRAN (INVOKED VIA MAKE GFORTRAN)
GNUFC = gfortran #You can get newer versions of gfortran, if you perform "scl enable devtoolset-10 bash" in your shell first
GNUQFC = /opt/rh/devtoolset-10/root/bin/gfortran
GNUFFLAGS = -O3 -ffast-math -march=native -p -opt-prefetch -fopenmp -std=legacy -llapack -cpp -J$(build) #Note that for new version of gfortran you might have to add -std=legacy or -fallow-argument-mismatch to compile random.f without errors!
#-fallow-argument-mismatch
GNUDBGFLAGS = -fcheck=bounds -fcheck=do -fcheck=mem -fcheck=pointer -p -O0 #-gdwarf-5 -O0 -Wall
#MPI VERSION (INVOKED VIA MAKE MPI)
MPIFC=mpif90
MPIFFLAGS = -fcx-limited-range -O3 -ffast-math -march=native -p -opt-prefetch -falign-commons -mcmodel=large -fopenmp -J$(build) -llapack -cpp -Dmpi_version #TODO: Check if all these flags are necessary!
#Syntax for running mpi calculations:
# - 1 machine with 12 cores: mpirun -np 12 genetic test.genetic
# - 4 machine with 12 cores: mpirun -np 48 --hostfile nodes.txt genetic test.genetic
# - nodes.txt specifies the nodes on which the program will run, the first mentioned note will perform the master thread
# - you have to start the calculation from the node with the master thread and have running sleep jobs for the other notes
# - TODO: Write a job file / submission script that automatizes this procedure
#mpirun -np 48 --hostfile nodes.txt genetic s_test-dist9-freeze.genetic > s_test-dist9-freeze.out &
######################################################################
#Extend search path for files (both .f and .incl files)
VPATH += $(src)
VPATH += $(src)parser
VPATH += $(src)parser/lib
VPATH += $(src)model
######################################################################
#Define objects for different Program parts (sorted in order of compilation)
parserlib_obj = strings.o long_keyread.o fileread.o keyread.o long_write.o
parser_obj = io_parameters.o accuracy_constants.o keys.o dim_parameter.o parameterkeys.o parse_errors.o parser.o
datamodule_obj = data_module.o #Compile this module before your model files and the genetic files
#model_obj = ptr_structure.o Potential_no3_5s_jcp2021_cart_corrected.o surface_mod.o matrix_form.o ctrans.o model.o weight.o adia.o
#model_obj = ptr_structure.o Potential_no3_5s_jcp2021_cart_corrected.o select_monom_mod.o maik_ctrans.o surface_mod.o matrix_form.o model.o weight.o adia.o
model_obj = ptr_structure.o ctrans.o surface_mod.o matrix_form.o model.o weight.o adia.o
mod_incl = mod_const.incl so_param.incl
random_obj = $(addprefix $(build), random.o)
genetic_obj = data_transform.o init.o write.o funcs.o marq.o lbfgsb.o idxsrt_mod.o fit_MeX.o mpi_fit_MeX.o genetic.o #content of data_transform and write is user specific, interfaces are fixed
objects = $(addprefix $(build), $(parserlib_obj) $(parser_obj) $(datamodule_obj) $(model_obj) $(genetic_obj) )
#plot_dip_obj = $(addprefix $(build), io_parameters.o accuracy_constants.o dim_parameter.oi model.o)
#Note: Since we are using modules, you have carefully choose the order of compilation and take dependencies between modules and subroutines into account!
######################################################################
# Lib path to pes libray
PATH_PES = $(HOME)/Documents/work/NO3/NO3_PES/NO3_PES_FABIAN/
PES_LIB = $(PATH_PES)libno3_pes_ffabian.a
FFLAGS += -I$(PATH_PES)
LDFLAGS = -L$(PATH_PES) -lno3_pes_ffabian
# define main goal
main = genetic
main1 = plot_dipole
.PHONY: ifort gfortran
ifort: $(main)
# define main compilation
gfortran: override FC = $(GNUFC)
gfortran: override FFLAGS = $(GNUFFLAGS)
gfortran: $(main)
$(main) : dirs $(random_obj) $(objects) $(PES_LIB)
$(FC) $(FFLAGS) $(random_obj) $(objects) $(LDFLAGS) -o $(bin)$(main)
parser.o : io_parameters.o keys.o dim_parameter.o parameterkeys.o parse_errors.o
$(FC) -c $(FFLAGS) $^ -o $@
$(build)%.o : %.f
$(FC) -c $(FFLAGS) $^ -o $@
$(build)%.o : %.f90
$(FC) -c $(FFLAGS) $^ -o $@
$(model_obj) : $(mod_incl)
######################################################################
# define name of additional recipes
.PHONY: clean neat remake debug test mpi gfortran gqfortran profile tar dirs
# define additionational recipes
trash= *__genmod* $(addprefix $(build),*__genmod* *.mod *.o)
clean:
$(RM) $(objects) $(trash)
neat: clean
$(RM) $(random_obj)
remake: clean $(main)
dirs:
@mkdir -p $(build) $(bin)
debug: override FFLAGS += $(DBGFLAGS)
debug: clean $(main)
cp $(infile) $(bin)
$(bin)$(main) $(bin)$(infile) | tee debug.out
modern: override FFLAGS = $(NEWFFLAGS)
modern: $(main)
gqfortran: override FC = $(GNUQFC)
gqfortran: override FFLAGS = $(GNUFFLAGS)
gqfortran: $(main)
gdebug: override FC = $(GNUFC)
gdebug: override FFLAGS = $(GNUFFLAGS) $(GNUDBGFLAGS)
gdebug: clean $(main)
mpi: override FC = $(MPIFC)
mpi: override FFLAGS = $(MPIFFLAGS)
mpi: $(main)
infile=hi-sing1-sig.genetic
gtest: override FC = $(GNUFC)
gtest: override FFLAGS = $(GNUFFLAGS)
gtest: clean $(main)
cp $(infile) $(bin)
$(bin)$(main) $(bin)$(infile) | tee test.out
gprofile: override FC = $(GNUFC)
gprofile: override FFLAGS = $(GNUFFLAGS) -pg
gprofile: clean $(main)
cp $(infile) $(bin)
test: clean $(main)
cp $(infile) $(bin)
$(bin)$(main) $(bin)$(infile) | tee test.out
profile: override FFLAGS += -pg
profile: clean test
date > profile
gprof $(bin)$(main) gmon.out >> profile
timestamp=$(shell date +"%FT%H-%M-%S")
tar:
date > INFO
tar --exclude-backups --exclude-vcs -czf tardir/geneticsrc_$(timestamp).tar src/ obj/ bin/ Makefile INFO

0
src/model/.giosaveVT1ebr Normal file
View File

38
src/model/JTmod.incl Normal file
View File

@ -0,0 +1,38 @@
!*** Relevant parameters for the analytic model
!*** offsets:
!*** offsets(1): morse equilibrium (N-H, Angström)
!*** offsets(2): reference angle (H-N-H)
!*** offsets(3): --
!*** pat_index: vector giving the position of the
!*** various coordinates (see below)
!*** ppars: polynomial parameters for tmcs
!*** vcfs: coefficients for V expressions.
!*** wzcfs: coefficients for W & Z expressions.
!*** ifc: inverse factorials.
integer matdim
parameter (matdim=5) ! matrix is (matdim)x(matdim)
real*8 offsets(2)
integer pat_index(maxnin)
! NH3 params
parameter (offsets=[2.344419d0,120.d0])
!##########################################################################
! coordinate order; the first #I number of coords are given to the
! ANN, where #I is the number of input neurons. The position i in
! pat_index corresponds to a coordinate, the value of pat_index(i)
! signifies its position.
!
! The vector is ordered as follows:
! a,xs,ys,xb,yb,b,rs**2,rb**2,b**2,
! es*eb, es**3, eb**3,es**2*eb, es*eb**2
! ri**2 := xi**2+yi**2 = ei**2; ei := (xi,yi), i = s,b
!
! parts not supposed to be read by ANN are marked by ';' for your
! convenience.
!##########################################################################
! a,rs**2,rb**2,es*eb,es**3,eb**3,es**2*eb,es*eb**2,b**2 #I=9 (6D)
parameter (pat_index=[1,2,3,4,5,6,7,8,9,10,11,12,13,14])
!**************************************************************************

120
src/model/adia.f90 Normal file
View File

@ -0,0 +1,120 @@
module adia_mod
implicit none
contains
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! % SUBROUTINE ADIA(N,P,NPAR,ymod,v,u,SKIP)
! %
! % determines the adiabatic energies by diagonalizing diabatic matrix.
! % The Eingenvalues are sorted according to the best fitting ordering
! % of the CI vectors.
! %
! % ATTENTION: The interface has changed. To sort by the ci's,
! % the datavalue of the current points are given
! %
! % input variables:
! % n: number of point (int)
! % p: parameter evector(double[npar])
! % npar: number of parameters (int)
! % skip: .false. if everything should be done
! %
! % output variables:
! % ymod: firtst nstat energies and than nci*ndiab ci's (double[ntot])
! % v: eigenvalues (double[ndiab])
! % u: eigenvectors (double[ndiab,ndiab])
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine adia(n,p,npar,ymod,vx,u,skip)
use dim_parameter,only: ndiab,nstat,ntot,nci,pst
use data_module,only: q_m,x1_m,x2_m,y_m
use diab_mod, only:diab
use data_matrix
!use dipole, only: diab
implicit none
integer i,j !running indices
integer iref !getting correction or refference
double precision e(ndiab,ndiab) !full diabatic matrix
double precision mx(ndiab,ndiab)
double precision my(ndiab,ndiab)
double precision vxs,vys,vxb,vyb
integer n !current point
integer npar !number of parameters
double precision p(npar) !parameters
double precision u(ndiab,ndiab),ut(ndiab,ndiab) !ci-vectors
double precision ymod(ntot) !fitted data
double precision vx(ndiab),vy(nstat) !eigen values
double precision,allocatable,dimension(:,:):: mat
logical skip,dbg
parameter (dbg=.false.)
double precision,dimension(2,2):: T,TT,TX,TY
! lapack variables
integer,parameter :: lwork = 1000
double precision work(lwork)
integer info
integer TYPES, BLK ! TYPE OF THE CALCULATION
! variabke for dgemm
double precision,dimension(ndiab,ndiab):: ex,ey,ez
double precision:: alpha
integer:: lda,ldb,beta,ldc
double precision,dimension(ndiab,ndiab):: temp1,temp2
call diab(ex,ey,ez,n,x1_m(:,n),x2_m(:,n),p)
! init eigenvector matrix
TYPES = int(p(pst(1,33)))
BLK = int(p(pst(1,33)+1)) ! BLOCK IF TYPE IS 3
u = 0.d0
vx=0.0d0
skip=.false.
ymod=0.0d0
if (TYPES .eq.1 ) then
! Trace of the potential
call trace_mat(ex,ey,ymod)
else if (TYPES .eq.2) then
! Eigenvalue decomposition of the potential
call Eigen(ex,ey,ymod)
else if (TYPES .eq.3) then
CALL BLOCK_DIAB(ex,ey,ymod,BLK)
else if (TYPES .EQ.4) then
call Full_diab_upper(ex,ey,ymod)
else if (TYPES .eq.5) then
call Transformation_mat(ex,vx,ymod)
ymod=0.0d0
else if (TYPES .eq.6) then
! transform the lz
call one_dia_upper(ez,ymod)
else
write(*,*) "Error in TYPE of calculation here",TYPES
stop
end if
if (dbg) then
do i=1,ndiab
write(*,'(5f14.6)') (ex(i,j),j=1,ndiab)
enddo
write(*,*)""
endif
end subroutine
subroutine matrix_mult(C,A,B,N)
implicit none
integer:: n,i,j,k
double precision,dimension(n,n):: A,B,C
do i = 1, n ! Rows of C
do j = 1, n ! Columns of C
C(i,j) = 0.0 ! Initialize element
do k = 1, n ! Dot product
C(i,j) = C(i,j) + A(i,k) * B(k,j)
end do
end do
end do
end subroutine
end module adia_mod

386
src/model/ctrans.f90 Normal file
View File

@ -0,0 +1,386 @@
module ctrans_mod
use dim_parameter, only: qn
contains
!! subroutine ctrans
subroutine ctrans(q,x1,x2)
implicit none
include 'nnparams.incl'
include 'JTmod.incl'
double precision,intent(in):: q(qn)
double precision,intent(out):: x1(qn),x2(qn)
double precision:: cart(3,4),qint(maxnin)
integer i
!cart(:,1)=0.0d0
!cart(1:3,2:4) = reshape([ q(4:12) ], shape(cart(1:3,2:4)))
cart(1,1)=q(1)
cart(2,1)=q(2)
cart(3,1)=q(3)
cart(1,2)=q(4)
cart(2,2)=q(5)
cart(3,2)=q(6)
cart(1,3)=q(7)
cart(2,3)=q(8)
cart(3,3)=q(9)
cart(1,4)=q(10)
cart(2,4)=q(11)
cart(3,4)=q(12)
call cart2int(cart,qint)
do i=1,qn
if (abs(qint(i)) .lt. 1.0d-5) qint(i) =0.0d0
enddo
x1(1:qn)=qint(1:qn)
!x1(2)=0.0d0
x1(5)=-x1(5)
x1(3)=-x1(3)
!x1(6)=0.0d0
x2(1:qn)=0.0d0 !qint(1:qn)
end subroutine ctrans
subroutine cart2int(cart,qint)
implicit none
! This version merges both coordinate transformation routines into
! one. JTmod's sscales(2:3) are ignored.
! This is the first version to be compatible with one of my proper 6D fits
! Time-stamp: <2024-10-22 13:52:59 dwilliams>
! Input (cartesian, in Angström)
! cart(:,1): N
! cart(:,1+i): Hi
! Output
! qint(i): order defined in JTmod.
! Internal Variables
! no(1:3): NO distances 1-3
! pat_in: temporary coordinates
! axis: main axis of NO3
include 'nnparams.incl'
include 'JTmod.incl'
real*8 cart(3,4),qint(maxnin)
real*8 no(3), r1, r2, r3
real*8 v1(3), v2(3), v3(3)
real*8 n1(3), n2(3), n3(3), tr(3)
real*8 ortho(3)
real*8 pat_in(maxnin)
logical ignore_umbrella,dbg_umbrella
logical dbg_distances
!.. Debugging parameters
!.. set umbrella to 0
parameter (ignore_umbrella=.false.)
! parameter (ignore_umbrella=.true.)
!.. break if umbrella is not 0
parameter (dbg_umbrella=.false.)
! parameter (dbg_umbrella=.true.)
!.. break for tiny distances
parameter (dbg_distances=.false.)
! parameter (dbg_distances=.true.)
integer k
!.. get N-O vectors and distances:
do k=1,3
v1(k)=cart(k,2)-cart(k,1)
v2(k)=cart(k,3)-cart(k,1)
v3(k)=cart(k,4)-cart(k,1)
enddo
no(1)=norm(v1,3)
no(2)=norm(v2,3)
no(3)=norm(v3,3)
!.. temporarily store displacements
do k=1,3
pat_in(k)=no(k)-offsets(1)
enddo
do k=1,3
v1(k)=v1(k)/no(1)
v2(k)=v2(k)/no(2)
v3(k)=v3(k)/no(3)
enddo
!.. compute three normal vectors for the ONO planes:
call xprod(n1,v1,v2)
call xprod(n2,v2,v3)
call xprod(n3,v3,v1)
do k=1,3
tr(k)=(n1(k)+n2(k)+n3(k))/3.d0
enddo
r1=norm(tr,3)
do k=1,3
tr(k)=tr(k)/r1
enddo
! rotate trisector
call rot_trisec(tr,v1,v2,v3)
!.. determine trisector angle:
if (ignore_umbrella) then
pat_in(7)=0.0d0
else
pat_in(7)=pi/2.0d0 - acos(scalar(v1,tr,3))
pat_in(7)=sign(pat_in(7),cart(1,2))
endif
!.. molecule now lies in yz plane, compute projected ONO angles:
v1(1)=0.d0
v2(1)=0.d0
v3(1)=0.d0
r1=norm(v1,3)
r2=norm(v2,3)
r3=norm(v3,3)
do k=2,3
v1(k)=v1(k)/r1
v2(k)=v2(k)/r2
v3(k)=v3(k)/r3
enddo
! make orthogonal vector to v3
ortho(1)=0.0d0
ortho(2)=v3(3)
ortho(3)=-v3(2)
!.. projected ONO angles in radians
pat_in(4)=get_ang(v2,v3,ortho)
pat_in(5)=get_ang(v1,v3,ortho)
pat_in(6)=dabs(pat_in(5)-pat_in(4))
!.. account for rotational order of atoms
if (pat_in(4).le.pat_in(5)) then
pat_in(5)=2*pi-pat_in(4)-pat_in(6)
else
pat_in(4)=2*pi-pat_in(5)-pat_in(6)
endif
pat_in(4)=rad2deg*pat_in(4)-offsets(2)
pat_in(5)=rad2deg*pat_in(5)-offsets(2)
pat_in(6)=rad2deg*pat_in(6)-offsets(2)
pat_in(7)=rad2deg*pat_in(7)
call genANN_ctrans(pat_in)
qint(:)=pat_in(:)
contains
!-------------------------------------------------------------------
! compute vector product n1 of vectors v1 x v2
subroutine xprod(n1,v1,v2)
implicit none
real*8 n1(3), v1(3), v2(3)
n1(1) = v1(2)*v2(3) - v1(3)*v2(2)
n1(2) = v1(3)*v2(1) - v1(1)*v2(3)
n1(3) = v1(1)*v2(2) - v1(2)*v2(1)
end subroutine
!-------------------------------------------------------------------
! compute scalar product of vectors v1 and v2:
real*8 function scalar(v1,v2,n)
implicit none
integer i, n
real*8 v1(*), v2(*)
scalar=0.d0
do i=1,n
scalar=scalar+v1(i)*v2(i)
enddo
end function
!-------------------------------------------------------------------
! compute norm of vector:
real*8 function norm(x,n)
implicit none
integer i, n
real*8 x(*)
norm=0.d0
do i=1,n
norm=norm+x(i)**2
enddo
norm=sqrt(norm)
end function
!-------------------------------------------------------------------
subroutine rot_trisec(tr,v1,v2,v3)
implicit none
real*8 tr(3),v1(3),v2(3),v3(3)
real*8 vrot(3)
real*8 rot_ax(3)
real*8 cos_phi,sin_phi
! evaluate cos(-phi) and sin(-phi), where phi is the angle between
! tr and (1,0,0)
cos_phi=tr(1)
sin_phi=dsqrt(tr(2)**2+tr(3)**2)
if (sin_phi.lt.1.0d-12) then
return
endif
! determine rotational axis
rot_ax(1) = 0.0d0
rot_ax(2) = tr(3)
rot_ax(3) = -tr(2)
! normalize
rot_ax=rot_ax/sin_phi
! now the rotation can be done using Rodrigues' rotation formula
! v'=v*cos(p) + (k x v)sin(p) + k (k*v) (1-cos(p))
! for v=tr k*v vanishes by construction:
! check that the rotation does what it should
call rodrigues(vrot,tr,rot_ax,cos_phi,sin_phi)
if (dsqrt(vrot(2)**2+vrot(3)**2).gt.1.0d-12) then
write(6,*) "ERROR: BROKEN TRISECTOR"
stop
endif
tr=vrot
call rodrigues(vrot,v1,rot_ax,cos_phi,sin_phi)
v1=vrot
call rodrigues(vrot,v2,rot_ax,cos_phi,sin_phi)
v2=vrot
call rodrigues(vrot,v3,rot_ax,cos_phi,sin_phi)
v3=vrot
end subroutine
!-------------------------------------------------------------------
subroutine rodrigues(vrot,v,axis,cos_phi,sin_phi)
implicit none
real*8 vrot(3),v(3),axis(3)
real*8 cos_phi,sin_phi
real*8 ortho(3)
call xprod(ortho,axis,v)
vrot = v*cos_phi + ortho*sin_phi+axis*scalar(axis,v,3)*(1-cos_phi)
end subroutine
!-------------------------------------------------------------------
real*8 function get_ang(v,xaxis,yaxis)
implicit none
! get normalized [0:2pi) angle from vectors in the yz plane
real*8 v(3),xaxis(3),yaxis(3)
real*8 phi
real*8 pi
parameter (pi=3.141592653589793d0)
phi=atan2(scalar(yaxis,v,3),scalar(xaxis,v,3))
if (phi.lt.0.0d0) then
phi=2*pi+phi
endif
get_ang=phi
end function
end subroutine cart2int
subroutine genANN_ctrans(pat_in)
implicit none
include 'nnparams.incl'
include 'JTmod.incl'
real*8 pat_in(maxnin)
real*8 raw_in(maxnin),off_in(maxnin),ptrans_in(7)
real*8 r0
real*8 a,b,xs,ys,xb,yb
integer k
off_in(1:7)=pat_in(1:7)
r0=offsets(1)
! transform primitives
! recover raw distances from offset coords
do k=1,3
raw_in(k)=off_in(k)+offsets(1)
enddo
do k=1,3
ptrans_in(k)=off_in(k)
enddo
! rescale ONO angles
ptrans_in(4)=deg2rad*off_in(4)
ptrans_in(5)=deg2rad*off_in(5)
ptrans_in(6)=deg2rad*off_in(6)
! rescale umbrella
ptrans_in(7)=off_in(7)*deg2rad
! compute symmetry coordinates
! A (breathing)
a=(ptrans_in(1)+ptrans_in(2)+ptrans_in(3))/dsqrt(3.0d0)
! ES
call prim2emode(ptrans_in(1:3),xs,ys)
! EB
call prim2emode(ptrans_in(4:6),xb,yb)
! B (umbrella)
b=ptrans_in(7)
! overwrite input with output
pat_in(pat_index(1))=a ! 1
pat_in(pat_index(2))=xs
pat_in(pat_index(3))=ys
pat_in(pat_index(4))=xb
pat_in(pat_index(5))=yb
pat_in(pat_index(6))=b
! totally symmetric monomials
pat_in(pat_index(7))=xs**2 + ys**2 ! 2
pat_in(pat_index(8))=xb**2 + yb**2 ! 3
pat_in(pat_index(9))=b**2 ! 9
pat_in(pat_index(10))=xs*xb+ys*yb ! 4
! S^3, B^3
pat_in(pat_index(11))=xs*(xs**2-3*ys**2) ! 5
pat_in(pat_index(12))=xb*(xb**2-3*yb**2) ! 6
! S^2 B, S B^2
pat_in(pat_index(13))=xb*(xs**2-ys**2) - 2*yb*xs*ys ! 7
pat_in(pat_index(14))=xs*(xb**2-yb**2) - 2*ys*xb*yb ! 8
do k=11,14
pat_in(pat_index(k))=tanh(0.1d0*pat_in(pat_index(k)))*10.0d0
enddo
end subroutine
subroutine prim2emode(prim,ex,ey)
implicit none
! Takes a 2D-vector prim and returns the degenerate modes x and y
! following our standard conventions.
real*8 prim(3),ex,ey
ex=(2.0d0*prim(1)-prim(2)-prim(3))/dsqrt(6.0d0)
ey=(prim(2)-prim(3))/dsqrt(2.0d0)
end
end module ctrans_mod

133
src/model/data_transform.f90 Normal file
View File

@ -0,0 +1,133 @@
! <subroutine for manipulating the input Data before the Fit
subroutine data_transform(q,x1,x2,y,wt,p,npar,p_act)
use accuracy_constants, only: dp,idp
use dim_parameter,only : nstat,pst,ntot,qn,numdatpt,ndiab,ndata,sets
use ctrans_mod, only: ctrans
use surface_mod, only: eval_surface
use data_matrix
! use david_ctrans_mod, only: ctrans_d
implicit none
! IN: variables
integer(idp),intent(in) :: npar
Real(dp),intent(in) :: q(qn,numdatpt)
Real(dp),intent(in) :: p(npar)
integer(idp),intent(in) :: p_act(npar)
! INOUT: variables
Real(dp),intent(inout) :: y(ntot,numdatpt)
Real(dp),intent(inout) :: wt(ntot,numdatpt)
! OUT: vairables
Real(dp), intent(out) :: x1(qn,numdatpt),x2(qn,numdatpt)
! internal variables
Real(dp),dimension(ndiab,ndiab)::mat_x,mat_y,mat_z,U,V
Real(dp),dimension(nstat) :: E
integer(idp) pt,i,j,k,l, TYPES, BLK ! types is for the type of calculation
! blk is for which block to fit
logical,parameter:: dbg = .false.
if (pst(2,33) .ne. 2) then
write(*,*) "Error in Paramater Keys, TYPE_CAL should be 2 parameter", pst(2,33)
stop
end if
TYPES = int(p(pst(1,33)))! TYPE OF THE CALCULATION
BLK= int(p(pst(1,33)+1))! BLOCK IF TYPE IS 3
write(*,*) "TYPE of calculation:",TYPES
pt=1
do i=1,sets ! loop over the number of sets
do j=1,ndata(i) ! loop over the nbr of points in each sets
! remember to increment pt at the end of the loop
call ctrans(q(1:qn,pt),x1(:,pt),x2(:,pt)) ! transform the coordinate
! get the reference U matrix
!if (j .eq. 3) then
! call eval_surface(E,V,U_ref,q(1:qn,pt))
! call transform_U(U_ref)
!endif
!do pt=1,numdatpt
!call ctrans(q(1:qn,pt),x1(:,pt),x2(:,pt))! ctrans the dipole cooordinate.
write(7,'(I3,*(E17.8))') pt,x1(:,pt)
call eval_surface(E,V,U,q(1:qn,pt))
! Transform U mmatrix
call transform_U(U) ! Transform the U matrix
! write U matrix on f16
if (dbg) then
!write(7,*) "U matrix at point", pt
do k=1,ndiab
write(50+i,'(2E17.8,5X,5E17.8)')x1(2:3,pt),(U(k,l),l=1,ndiab)
enddo
write(50+i,*) ""
endif
!call overlap(U_ref,U)
call Y2mat(y(1:ntot,pt),mat_x,mat_y,mat_z)
if (TYPES .eq.1 ) then
! Trace of the potential
call trace_mat(mat_x,mat_y,y(1:ntot,pt))
else if (TYPES .eq.2) then
! Eigenvalue decomposition of the potential
call Eigen(mat_x,mat_y,y(1:ntot,pt))
else if (TYPES .eq.3) then
! Adiabatic transformation of the potential
call adiabatic_transform(mat_x,U)
call adiabatic_transform(mat_y,U)
call block_diab(mat_x,mat_y,y(1:ntot,pt),BLK)
else if (TYPES .eq.4) then
! Write the full upper diabatic matrix
call adiabatic_transform(mat_x,U)
call adiabatic_transform(mat_y,U)
! and write the full diabatic matrix to y
! This is the full diabatic matrix
call Full_diab_upper(mat_x,mat_y,y(1:ntot,pt))
else if (TYPES .eq.5) then
!call adiabatic_transform(mat_x,U)
!call adiabatic_transform(mat_y,U)
call Transformation_mat(U,E,y(1:ntot,pt))
else if (TYPES .eq.6) then
! Just do the adiabatic transformation and write the matrix
! transform the lz
call adiabatic_transform(mat_z,U)
call one_dia_upper(mat_z,y(1:ntot,pt))
else
write(*,*) "Error in TYPE of calculationss",TYPES
write(*,*) "the value:,", p(pst(1,33))
stop
end if
pt=pt+1
enddo ! j
write(34,*) "#---- End of set ", i
write(7,*) "#---- End of set ", i
enddo ! i
!enddo
call weight(wt,y)
end subroutine

84
src/model/keys.f90 Normal file
View File

@ -0,0 +1,84 @@
module keys_mod
implicit none
contains
!program gen_key
! implicit none
! call init_keys()
!end program gen_key
subroutine init_keys
use io_parameters, only: key
character(len=1) prefix(4)
parameter (prefix=['N','P','A','S'])
!character (len=20) key(4,25)
integer,parameter:: np=33
character(len=16) parname(np)
integer i,j
! Defining keys for potential
! the electronic state of NO3 A2' E" E'
! Naming convention
! the keys for Lx and Ly
! the coupling between A2' and A2"
parname(1)='LXYVA2O1'
parname(2)='LXYVE1O1'
parname(3)='LXYVE2O1'
parname(4)='LXYVA2O2'
parname(5)='LXYVE1O2'
parname(6)='LXYVE2O2'
! W & Z of E1
parname(7)='LXYWZE1O0'
parname(8)='LXYWZE1O1'
parname(9)='LXYWZE1O2'
parname(10)='LXYWZE2O0'
parname(11)='LXYWZE2O1'
parname(12)='LXYWZE2O2'
! WW and Z Pseudo between E1 and E2
! p STANDS FOR PSEUDO JAHN-TELLER
parname(13)='LXYPE1E2O0'
parname(14)='LXYPE1E2O1'
parname(15)='LXYPE1E2O2'
! no order 3
! PSEUDO A2 & E1
parname(16)='LXYPA2E1O0'
parname(17)='LXYPA2E1O1'
parname(18)='LXYPA2E1O2'
! Pseudo JAHN-TELLER BETWEEN A2 AND E1
parname(19)='LXYPA2E2O0'
parname(20)='LXYPA2E2O1'
parname(21)='LXYPA2E2O2'
! keys for lz
parname(22)='LZWZE1O1'
parname(23)='LZWZE1O2'
parname(24)='LZWZE2O1'
parname(25)='LZWZE2O2'
parname(26)='LZPE1E2O0'
parname(27)='LZPE1E2O1'
parname(28)='LZPE1E2O2'
parname(29)='LZPA2E1O1'
parname(30)='LZPA2E1O2'
parname(31)='LZPA2E2O1'
parname(32)='LZPA2E2O2'
parname(33)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE
do i=1,np
do j=1,4
key(j, i)=prefix(j)//trim(parname(i))//':'
write(8,*) key(j,i)
enddo
write(8,*) ''
enddo
end subroutine
end module keys_mod

374
src/model/matrix_form.f90 Normal file
View File

@ -0,0 +1,374 @@
module data_matrix
use dim_parameter, only:ndiab,nstat,ntot,pst
! use surface_mod, only: eval_surface
contains
! subroutine trace
subroutine trace_mat(mx,my,y)
IMPLICIT NONE
integer::i
double precision,intent(inout):: y(:)
double precision, intent(in):: mx(:,:),my(:,:)
y=0.0d0
!y(1)=mx(4,4)+mx(5,5)
do i=1,ndiab
y(1)=y(1)+mx(i,i)
y(2)=y(2)+my(i,i)
enddo
END SUBROUTINE trace_mat
!! subroutine Ydata to matrix
subroutine Y2mat(Y,Mx,My,mz)
IMPLICIT NONE
integer:: ii,i,j
double precision, intent(in):: y(:)
double precision,dimension(ndiab,ndiab),intent(out):: mx, my,mz
!if (ndiab .ne. 4 ) then
!write(*,*) " NDIAB should be equal to 4",NDIAB
!write(*,*) "CHECK DATA_TRANSFORM TO MAKE IT ADAPTABLE"
!stop
!endif
ii=1
do i=1,ndiab
do j=1,i
! !mx
mx(i,j)=y(ii)
! ! My
my(i,j)=y( (ntot/3)+ii)
! remember to adjust here I added the energy
mz(i,j)= y(2*(ntot/3)+ ii )
!
ii=ii+1
enddo
enddo
call copy_2_upper(mx)
call copy_2_upper(my)
call copy_2_upper(mz)
end subroutine
subroutine Full_diab_upper(mx,my,y)
implicit none
double precision,intent(inout) :: y(:)
double precision, intent(in) :: mx(ndiab,ndiab), my(ndiab,ndiab)
integer i,j,ii
ii=1
y=0.0d0
do i=1,ndiab
do j=i,ndiab
! mx
y(ii) = mx(i,j)
! my
y((ntot/2)+ii) = my(i,j)
! increment the index
ii=ii+1
enddo
enddo
end subroutine Full_diab_upper
subroutine one_dia_upper(m,y)
implicit none
double precision,intent(inout) :: y(:)
double precision, intent(in) :: m(ndiab,ndiab)
integer i,j,ii
ii=1
y=0.0d0
do i=1,ndiab
do j=i,ndiab
! mx
y(ii) = m(i,j)
! increment the index
ii=ii+1
enddo
enddo
end subroutine one_dia_upper
Subroutine adiabatic_transform(mx,U)
implicit none
double precision, intent(inout) :: mx(ndiab,ndiab)
double precision, dimension(:,:), intent(inout) :: U
double precision, dimension(ndiab,ndiab) :: temp1, temp2
integer i, j
!call transform_U(U) ! Transform the U matrix
! Transform mx and my to adiabatic basis
temp1 = matmul(mx, transpose(U))
mx = matmul(U, temp1)
!temp2 = matmul(my, transpose(U))
!my = matmul(U, temp2)
end subroutine adiabatic_transform
! the eigenvalue of the dipole
SUBROUTINE Eigen(mx,my,Yres)
implicit none
double precision,dimension(:,:),intent(inout) :: mx,my
double precision,dimension(:),intent(out) :: Yres
double precision,dimension(ndiab) :: vx,vy
double precision,dimension(size(mx,1),size(my,2)) :: temp
! create a temorary matrix fo the eigenvctors
double precision, allocatable :: mux(:,:), muy(:,:)
! Lapak parameters
integer :: n,info,i
integer,parameter :: lwork = 100
double precision :: work(lwork)
! temporary
double precision:: max_row
Yres = 0.0d0
Allocate(mux,source=mx)
call DSYEV('V', 'U', size(mx,1), mux, size(mx,1), vx, work, lwork, info)
mx=mux
if (info /= 0) then
write(*,*) "Error in Eigenvalue decomposition of mx info = ", info
stop
end if
deallocate(mux)
Allocate(muy,source=my)
call DSYEV('V', 'U', size(my,1), muy, size(my,1), vy, work, lwork, info)
if (info /= 0) then
write(*,*) "Error in Eigenvalue decomposition of my info = ", info
stop
end if
deallocate(muy)
Yres(1:size(mx,1)) = vx(1:size(mx,1))
do i=1,size(mx,1)
max_row=maxloc(abs(mx(:,i)),1)
!yres(size(mx,1)+i)=(mx(max_row,i))**2
!yres(size(mx,1)+i)=real(max_row)
enddo
!Yres(size(mx,1)+1:2*size(mx,1)) = vy(1:size(my,1))
end subroutine
subroutine copy_2_upper(m)
implicit none
double precision, intent(inout) :: m(:,:)
integer :: i,j
! copy the lower part of the matrix to the upper part
do i=1,size(m,1)
do j=1,i-1
m(j,i) = m(i,j)
enddo
enddo
end subroutine copy_2_upper
subroutine coppy_2_low(m)
implicit none
double precision, intent(inout) :: m(:,:)
integer :: i,j
! copy the upper part of the matrix to the lower part
do i=1,size(m,1)
do j=i+1,size(m,2)
m(j,i) = m(i,j)
enddo
enddo
end subroutine coppy_2_low
!1 SUBROUTNE BLOCKS
!! EACH BLOCK OF dIABTIC MATRIX
SUBROUTINE block_diab(mx,my,Y,block)
implicit none
double precision, intent(inout):: Y(:)
double precision, intent(in) :: mx(ndiab,ndiab), my(ndiab,ndiab)
integer, intent(in) :: block
integer i,j,ii,nn
y=0.0d0
select case (block)
case(1)
! fill the first E1 block state 2 &3
y(1)=mx(2,2)
y(2)=mx(2,3)
!y(3)=mx(3,2)
y(4)=mx(3,3)
!y(5)=my(2,2)
!y(6)=my(2,3)
!y(7)=my(3,2)
!y(8)=my(3,3)
case(2)
! fill the second E2 block state 4 & 5
y(1)=mx(4,4)
y(2)=mx(4,5)
!y(3)=mx(5,4)
y(4)=mx(5,5)
y(5)=my(4,4)
y(6)=my(4,5)
!y(7)=my(5,4)
y(8)=my(5,5)
case(3)
! Filling the pseudo block E1 and E2
y(1)=mx(2,4)
y(2)=mx(2,5)
y(3)=mx(3,4)
y(4)=mx(3,5)
y(5)=my(2,4)
y(6)=my(2,5)
y(7)=my(3,4)
y(8)=my(3,5)
case(4)
! filling the block of A2 coupling with E1
y(1)=mx(1,2)
y(2)=mx(1,3)
y(3)=mx(2,1)
y(4)=mx(3,1)
!y(5)=my(1,2)
!y(6)=my(1,3)
!y(7)=my(2,1)
!y(8)=my(3,1)
case(5)
! couplinng A2 with E2
Y(1)=mx(1,4)
Y(2)=mx(1,5)
!Y(3)=mx(4,1)
!Y(4)=mx(5,1)
Y(5)=my(1,4)
Y(6)=my(1,5)
!Y(7)=my(4,1)
!Y(8)=my(5,1)
case(6)
! Filling A only
y(1)=mx(1,1)
y(5)=my(1,1)
case default
write(*,*) "Error in block_diab subroutine, block not recognized"
write(*,*) "The block is:", block
stop
end select
end subroutine block_diab
subroutine ident(A)
implicit none
integer i,j
double precision,intent(inout)::A(:,:)
do i=1,size(A,1)
do j=1,size(A,1)
if (i==j) then
A(i,j)=1.0d0
else
A(i,j)=0.0d0
endif
enddo
enddo
end subroutine
! subroutine trasform the U matrix
subroutine transform_U(U)
implicit none
double precision, intent(inout) :: U(ndiab,ndiab)
double precision :: U_ref(ndiab,ndiab), V(ndiab,ndiab), E(nstat)
integer i,max_row
double precision:: dot_prod,q_ref(9)
logical,parameter:: dbg_sign =.true.
!q_ref= [1.000174,0.000000,0.000000,-0.503595,-0.872253,0.000000,-0.530624,0.919068,0.000000]
!call eval_surface(E,V,U_ref,q_ref,p) ! get the reference transformation matrix
do i=1,ndiab
max_row = maxloc(abs(U(:,i)),1)
if (U(max_row,i) .lt. 0) then
U(:,i) = -1*U(:,i)
endif
enddo
!dot_prod=dot_product(U(2:3,4),U_ref(2:3,4))
!if (dot_prod .lt. 0.0d0) then
! U(:,4) = -1.0d0*U(:,4)
!endif
end subroutine transform_U
subroutine write_type_calc(p,id_write)
! Subroutine to write the type of calculation
implicit none
double precision, intent(in) :: p(:)
integer, intent(in) :: id_write
integer :: type_calc, blk
type_calc = int(p(pst(1,33)))
blk = int(p(pst(1,33)+1))
if (type_calc ==1) then
write(id_write,*) "Type of calculation: TRACE"
else if (type_calc ==2) then
write(id_write,*) "Type of calculation: EIGENVALUE"
else if (type_calc ==3) then
IF (blk == 1) then
write(id_write,*) "Type of calculation: E1 BLOCK"
ELSE IF (BLK ==2) THEN
write(id_write,*) "Type of calculation: E2 BLOCK"
ELSE IF (BLK ==3) THEN
write(id_write,*) "Type of calculation: Pseudo E1 and E2 BLOCK"
ELSE IF (BLK ==4) THEN
write(id_write,*) "Type of calculation: COUPLING A2 with E1 BLOCK"
ELSE IF (BLK ==5) THEN
write(id_write,*) "Type of calculation: COUPLING A2 with E2 BLOCK"
ELSE IF (BLK ==6) THEN
write(id_write,*) "Type of calculation: A2 ONLY"
ELSE
write(id_write,*) "Type of calculation: Diabatic transformation with unknown block size", blk
END IF
else if (type_calc ==4) then
write(id_write,*) "Type of calculation: Full Diabatic Matrix"
else if (type_calc ==5) then
write(id_write,*) "Type of calculation: Transformation matrix U"
else
write(id_write,*) "Error in type of calculation:", type_calc
stop
end if
END SUBROUTINE write_type_calc
!! subroutine for writting the transformtion matrix U
subroutine Transformation_mat(temp,v,y)
implicit none
double precision, intent(in) :: temp(ndiab,ndiab), v(:)
double precision, intent(inout) :: y(:)
double precision :: U(ndiab,ndiab )
integer i,j,ii
U(1:ndiab,1:ndiab) = temp(1:ndiab,1:ndiab)
!call transform_U(U,P)
y=0.0d0
!y(1:4) = v(1:4) ! copy the first 4 elements of v to y
ii=1
do i=1,ndiab
do j=1,ndiab
y(ii) = U(i,j)
ii=ii+1
enddo
enddo
y(ii:30)=v(:)
end subroutine
! compute the overlap between U matrix
subroutine overlap(U_ref,U)
implicit none
double precision, intent(in):: U_ref(ndiab,ndiab)
double precision, intent(inout):: U(ndiab,ndiab)
double precision:: over
integer i
do i=1,ndiab
over=dot_product(U_ref(:,i),U(:,i))
if (over .lt. 0.0d0 ) then
U(:,i)=-U(:,i)
endif
enddo
end subroutine
end module

507
src/model/model.f90 Normal file
View File

@ -0,0 +1,507 @@
! Author: jnshuti
! Created: 2025-10-03 14:09:49
! Last modified: 2025-10-03 14:10:10 jnshuti
! model for L-matrix of NO3 radical
module diab_mod
use accuracy_constants, only: dp, idp
use dim_parameter, only: ndiab, nstat, ntot,qn,pst
implicit none
private
public :: diab
contains
subroutine diab(lx,ly,lz,n,x1,x2,p)
implicit none
real(dp), intent(out),dimension(ndiab,ndiab):: lx,ly,lz
real(dp), intent(in), dimension(qn):: x1,x2
real(dp), intent(in),dimension(:):: p
integer(idp),intent(in):: n
call Lx_diab(lx,x1,x2,p)
call Ly_diab(ly,x1,x2,p)
call Lz_diab(lz,x1,x2,p)
end subroutine diab
subroutine Lx_diab(E,q,t,p)
implicit none
real(dp),dimension(ndiab,ndiab), intent(out):: E
real(dp),dimension(:),intent(in):: q,t
real(dp),dimension(:),intent(in):: p
real(dp):: xs,ys,xb,yb,a,b
real(dp):: v3_vec(8), v2(6)
integer(idp):: i,j,id
! check the dimension of the matrix
if (size(E,1) .ne. ndiab) then
write(*,*) " Error in Lx_diab: wrong dimension of L matrix ", size(E,1)
stop
endif
! rewrite the coordinate array q into symmetry adapted coordinates
call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
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
e = 0.0_dp
id = 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) & ! 3 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)
! W and Z term of E1
! order 0
id=id+1 ! 7
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 ! 8 ! 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 ! 9 ! 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
! 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 ! 10
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 !112 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 ! 12 ! 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
! make the dipole E = b* E
e = b * e
! E1 X E2
! WW and ZZ
id =id+1 ! 13
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 ! 14 ! 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 ! 15
do i=1,3 ! 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 ! 16
e(1,3)=e(1,3)+b*(p(pst(1,id)))
! order 1
id = id +1 ! 17
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 ! 18
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 !19
e(1,5)=e(1,5)+p(pst(1,id))
! order 1
id = id +1 ! 20
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 ! 21
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)
e(1,4:5) = b* e(1,4:5)
call copy_2_lower_triangle(e)
end subroutine Lx_diab
! Ly matrix
subroutine Ly_diab(e,q,t,p)
implicit none
real(dp),dimension(ndiab,ndiab), intent(out):: e
real(dp),dimension(:),intent(in):: q,t
real(dp),dimension(:),intent(in):: p
real(dp):: xs,ys,xb,yb,a,b
real(dp):: v2(6)
integer(idp):: i,j,id
! check the dimension of the matrix
if (size(e,1) .ne. ndiab) then
write(*,*) " Error in Ly_diab: wrong dimension of L matrix ", size(e,1)
stop
endif
! rewrite the coordinate array q into symmetry adapted coordinates
call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
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
e = 0.0_dp
! 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))
! W and Z of E1
! order 0
id=id+1 ! 7
e(2,3)=e(2,3)+p(pst(1,id))
! order 1
id=id+1 ! 8
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 ! 9
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
!! W and Z of E2
!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! order 0
id=id+1 ! 10
e(4,5)=e(4,5)+p(pst(1,id))
! order 1
id=id+1 ! 11
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 ! 12
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
! PSEUDO JAHN-TELLER E1 AND E2
e = b* e
!ORDER 0
id=id+1 ! 13
e(2,5)=e(2,5)+p(pst(1,id))
e(3,4)=e(3,4)+p(pst(1,id))
! order 1
id=id+1 ! 14
e(2,4)=e(2,4)+((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)+((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)+((-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)+((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 ! 15
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)
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)
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)
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)
! no order 3
!!!!!!!!!!!!!!!!
! the coupling A2 & E1
! #####################
! order 0
id=id+1 ! 16
e(1,2)=e(1,2)+(p(pst(1,id)))
! order 1
id=id+1 ! 17
e(1,2)=e(1,2)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
e(1,3)=e(1,3)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
! order 2
id=id+1 !18
e(1,2)=e(1,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(1,3)=e(1,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))
! COUPLING OF A2 WITH E2
!#######################################################################################
!###############################################################################
! order 0
id = id+1 !19
e(1,4)=e(1,4)+p(pst(1,id))
! order 1
id=id+1 ! 20
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 ! 21
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
e(1:4,5) = b * e(1:4,5)
call copy_2_lower_triangle(e)
end subroutine Ly_diab
! Lz matrix
subroutine Lz_diab(e,q,t,p)
implicit none
real(dp),dimension(ndiab,ndiab), intent(out):: e
real(dp),dimension(:),intent(in):: q,t
real(dp),dimension(:),intent(in):: p
real(dp):: xs,ys,xb,yb,a,b
real(dp):: v2(6)
integer(idp):: i,j,id
! check the dimension of the matrix
if (size(e,1) .ne. ndiab) then
write(*,*) " Error in Lz_diab: wrong dimension of e matrix ", size(e,1)
stop
endif
call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
e = 0.0_dp
! id for lz
id = 22 ! has to be
! the diagonal terms
! the v-term is 0th order and 3rd order.
! There is no zeroth order for diagonal
! w and z of E''
! order 1
id = id ! 22
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 ! 23
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
! W and Z of E'
! order 1
id = id +1 ! 24
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 ! 25
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
! the coupling
! Pseudo of E' and E''
! it must have odd power of b
id = id +1 !26
! order 0
e(2,4) = e(2,4)
e(3,5) = e(3,5)
e(2,5) = e(2,5) + b*(p(pst(1,id)))
e(3,4) = e(3,4) - b*(p(pst(1,id)))
! order 1
id = id +1 !27
e(2,4) = e(2,4) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
e(3,5) = e(3,5) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
e(2,5) = e(2,5) - b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
e(3,4) = e(3,4) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
! order 2
id = id +1 !28
do i=1,3
e(2,4) = e(2,4) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
e(3,5) = e(3,5) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
e(2,5) = e(2,5) + b*(p(pst(1,id)+(i-1)))*v2(i)
e(3,4) = e(3,4) - b*(p(pst(1,id)+(i-1)))*v2(i)
enddo
! no third order
! the coupling between A2' and E''
! order 1
id = id +1 !29
e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)*xb))
e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)*yb))
id = id +1 !30
! order 2
do i=1,3
e(1,2) = e(1,2) + b*(p(pst(1,id)+(i-1)))*v2(i)
e(1,3) = e(1,3) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
enddo
! the coupling of A2' and E'
! order 1
id = id +1 !31
e(1,2) = e(1,2) + (p(pst(1,id))*xs + p(pst(1,id)*xb))
e(1,3) = e(1,3) - (p(pst(1,id))*ys + p(pst(1,id)*yb))
id = id +1 ! 32
! order 2
do i=1,3
e(1,2) = e(1,2) + (p(pst(1,id)+(i-1)))*v2(i)
e(1,3) = e(1,3) + (p(pst(1,id)+(i-1)))*v2(i+3)
enddo
call copy_2_lower_triangle(e)
end subroutine Lz_diab
subroutine rewrite_coord(q,a,xs,ys,xb,yb,b,start)
implicit none
real(dp),dimension(:),intent(in):: q
real(dp),intent(out):: xs,ys,xb,yb,a,b
integer(idp),intent(in):: start
integer(idp):: i,j
a= q(start)
xs = q(start+1)
ys = q(start+2)
xb = q(start+3)
yb = q(start+4)
b = q(start+5)
end subroutine rewrite_coord
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 diab_mod

43
src/model/nnparams.incl Normal file
View File

@ -0,0 +1,43 @@
!**** Declarations
real*8 pi
real*8 hart2eV, eV2hart
real*8 hart2icm, icm2hart
real*8 eV2icm, icm2eV
real*8 deg2rad, rad2deg
integer maxnin,maxnout
!**********************************************************
!**** Parameters
!*** maxnin: max. number of neurons in input layer
!*** maxnout: max. number of neurons in output layer
parameter (maxnin=14,maxnout=15)
!**********************************************************
!**** Numerical Parameters
!*** infty: largest possible double precision real value.
!*** iinfty: largest possible integer value.
! 3.14159265358979323846264338327950...
parameter (pi=3.1415926536D0)
!**********************************************************
!**** Unit Conversion Parameters
!*** X2Y: convert from X to Y.
!***
!*** hart: hartree
!*** eV: electron volt
!*** icm: inverse centimeters (h*c/cm)
!****
!*** deg: degree
!*** rad: radians
parameter (hart2icm=219474.69d0)
parameter (hart2eV=27.211385d0)
parameter (eV2icm=hart2icm/hart2eV)
parameter (icm2hart=1.0d0/hart2icm)
parameter (eV2hart=1.0d0/hart2eV)
parameter (icm2eV=1.0d0/eV2icm)
parameter (deg2rad=pi/180.0d0)
parameter (rad2deg=1.0d0/deg2rad)

85
src/model/surface_mod.f90 Normal file
View File

@ -0,0 +1,85 @@
module surface_mod
use accuracy_constants, only: dp
implicit none
private
public eval_surface
contains
subroutine eval_surface(e, w, u, x1)
use accuracy_constants, only: dp, idp
use dim_parameter, only: ndiab
implicit none
real(dp), dimension(:, :), intent(out) :: w, u
real(dp), dimension(:), intent(out) :: e
real(dp), dimension(:), intent(in) :: x1
real(dp), allocatable, dimension(:, :) :: Mat
! debug parameter
logical, parameter:: dbg=.false.
integer(kind=idp):: i,j
! lapack variables
integer(kind=idp), parameter :: lwork = 1000
real(kind=dp) work(lwork)
integer(kind=idp) info
!write(*,*)"# Calling the potential routine "
call init_pot_para
call potentialno35s(W,X1)
allocate (Mat, source=w)
call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info)
if( info .ne. 0) then
write(*,*) " Error in eigenvalues decomposition routine of potential info=", info
stop
endif
u(:, :) = Mat(:, :)
deallocate (Mat)
if (dbg) then
do i=1,ndiab
write(19,99) e(i),(U(i,j),j=1,ndiab)
enddo
write(19,*)""
endif
99 format(2x,f16.8,2X,5f16.8)
end subroutine eval_surface
! subroutine init_surface(p)
! use dim_parameter, only: ndiab, nstat, ntot, nci ,qn
! use parameterkeys, only: parameterkey_read
! use fileread_mod, only: get_datfile, internalize_datfile
! use io_parameters, only: llen
! use accuracy_constants, only: dp
! implicit none
! real(dp), dimension(:), allocatable, intent(out) :: p
! character(len=llen), allocatable, dimension(:) :: infile
!
! qn = 9
! ndiab = 4
! nstat = 4
! nci = 4
! ntot = ndiab + nstat + nci
!
! block
! character(len=:),allocatable :: datnam
! integer :: linenum
! !get parameter file
! call get_datfile(datnam)
! !internalize datfile
! call internalize_datfile(datnam, infile, linenum, llen)
! end block
!
! !read parameters from file
! block
! real(dp), dimension(:), allocatable :: p_spread
! integer,dimension(:),allocatable :: p_act
! integer :: npar
! real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp
! call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread)
! end block
! end subroutine init_surface
end module surface_mod

50
src/model/weight.f Normal file
View File

@ -0,0 +1,50 @@
! <Subroutine weight(wt,y,ntot,numdatpt)
subroutine weight(wt,y)
use dim_parameter, only: nstat,ndiab,nci,ntot,numdatpt,
> hybrid,wt_en2ci,wt_en,wt_ci
implicit none
! data arrays and their dimensions
double precision wt(ntot,numdatpt),y(ntot,numdatpt)
! loop index
integer i,j,k,n
do i=1,numdatpt
wt(1,i)=1.d0
enddo
call norm_weight(wt,ntot,numdatpt)
end
!----------------------------------------------------------------------------------------------------
! <Subroutine norm_weight(wt,ntot,numdatpt)
subroutine norm_weight(wt,ntot,numdatpt)
implicit none
integer ntot,numdatpt
double precision norm,wt(ntot,numdatpt)
integer i,j,count
write(6,*) 'Normalizing Weights...'
norm=0.d0
count = 0
do i=1,numdatpt
do j=1,ntot
norm = norm + wt(j,i)*wt(j,i)
if (wt(j,i).gt.0.d0) count=count+1
enddo
enddo
norm = dsqrt(norm)
if(norm.gt.0.d0) then
do i=1,numdatpt
do j=1,ntot
wt(j,i) = wt(j,i)/norm
enddo
enddo
else
write(6,*) 'Warning: Norm of Weights is Zero'
endif
Write(6,'(''No. of weigthed data points:'',i0)') count
end subroutine

763
src/model/write.f Normal file
View File

@ -0,0 +1,763 @@
module write_mod
implicit none
! unit conversion
double precision ,parameter :: h2icm = 219474.69d0
double precision, parameter :: au2Debye = 2.541746d0
character(len=250), parameter :: sep_line = '(250("-"))'
character(len=250), parameter :: block_line = '(250("="))'
contains
! <Subroutine for writing the Output
subroutine write_output
> (q,x1,x2,y,wt,par,p_act,p_spread,nset,npar,
> flag,lauf)
use adia_mod, only: adia
use dim_parameter,only: qn,ntot,numdatpt,ndiab
use ctrans_mod,only: ctrans
implicit none
! IN: variables
integer lauf
integer flag !< 0= initial output 1=fit not converged 2= Fit Converged, 3= max iteration reached
integer npar,nset
double precision par(npar,nset),p_spread(npar)
integer p_act(npar)
double precision q(qn,numdatpt),x1(qn,numdatpt),x2(qn,numdatpt)
double precision y(ntot,numdatpt),wt(ntot,numdatpt)
! INTERNAL: Variables
integer,parameter :: id_out = 20 , std_out = 6
integer pt
integer i, id_print
double precision, allocatable :: ymod(:,:)
double precision, allocatable :: ew(:,:)
double precision, allocatable :: ev(:,:,:)
logical skip
allocate(ymod(ntot,numdatpt))
allocate(ew(ndiab,numdatpt))
allocate(ev(ndiab,ndiab,numdatpt))
skip=.false.
! get Model Outputs for all geometries for current best parameter set par(:,1)
do pt=1,numdatpt
call adia(pt,par(1:npar,1),npar,ymod(1:ntot,pt),
> ew(1:ndiab,pt),ev(1:ndiab,1:ndiab,pt),skip)
call ctrans(q(:,pt),x1(:,pt),x2(:,pt))
enddo
! Initial write print everything you want to see before the fit and return
if(flag.eq.0) then
call print_parameterstate(std_out,par(:,1),p_act,npar)
call print_ErrorSummary(std_out,y,ymod,wt)
! print Data into the plotfiles
return
endif
! open output files for individual makro iterations
call open_outfile(id_out,lauf)
! print Data into the plotfiles
call print_plotfiles(x1,y,wt,ymod)
! print Genetic output into files
do i=1, 2
if (i.eq.1) then
id_print= std_out
else
id_print= id_out
endif
write(id_print,'("Writing Iteration: ",i4)') lauf
write(id_print,block_line)
! write data information only in outfile
if(i.eq.2) then
call print_data(id_print,x1,y,ymod,wt)
call print_Set_Errors(id_print,y,ymod,wt)
endif
call print_parameterblock
> (id_print,par(:,1),p_act,p_spread,npar)
call print_ErrorSummary(id_print,y,ymod,wt)
enddo
call print_fortranfile(par(:,1),npar)
! write the type of calc at the end of the output
close (id_out)
deallocate(ymod,ev,ew)
end subroutine
!----------------------------------------------------------------------------------------------------
! <subroutine for scan seperated Error analysis>
subroutine print_Set_Errors(id_out,y, ymod, wt)
use io_parameters,only: llen
use dim_parameter,only: ndata,nstat,ntot,numdatpt,sets
integer , intent(in) :: id_out
double precision, intent(in) :: y(ntot,numdatpt),
> ymod(ntot,numdatpt), wt(ntot,numdatpt)
integer :: set, setpoint, pt
double precision :: Set_rms(sets,ntot), Set_num(sets,ntot)
double precision :: Total_rms, Total_Energy_rms,Energy_rms(nstat)
character(len=llen) fmt
write(id_out,'(A)') 'Errors in icm for individual Sets' //
> '(specified by sets: and npoints:)'
write(id_out,'(A5,3A16)')'Set','Total',
> 'Total_Energy', 'Energy[nstat]'
write(id_out,sep_line)
write(fmt,'("(I5,2f16.1,",I2,"f16.1)")') nstat
Set_rms = 0.d0
pt = 0
do set=1, sets
do setpoint=1, ndata(set)
pt = pt + 1
where(wt(:,pt) > 0.d0)
Set_rms(set,:) = Set_rms(set,:)+(ymod(:,pt)-y(:,pt))**2
Set_num(set,:) = Set_num(set,:) + 1
end where
enddo
Total_rms
> = dsqrt(sum(Set_rms(set,:))
> / (sum(Set_num(set,:))))
Total_Energy_rms
> = dsqrt(sum(Set_rms(set,1:nstat))
> / (sum(Set_num(set,1:nstat))))
Energy_rms(1:nstat)
> = dsqrt(Set_rms(set,1:nstat)
> / (Set_num(set,1:nstat)))
write(id_out,fmt) set, Total_rms*h2icm, Total_Energy_rms*h2icm,
> Energy_rms(1:nstat)*h2icm
enddo
write(id_out,block_line)
write(id_out,*) ''
end subroutine print_Set_Errors
!----------------------------------------------------------------------------------------------------
! <subroutine for printing the parameter and the pst vector in fortran readable style for including the fitted parameters in other programs
subroutine print_fortranfile(p,npar)
use io_parameters,only: maxpar_keys
use dim_parameter,only: pst
implicit none
! IN: variables
integer npar
double precision p(npar)
! INTERNAL: variables
integer i
integer, parameter :: id_out = 49
character(len=32), parameter :: fname ='fit_genric_bend_no3.f90'
open(id_out,file=fname)
30 format(6x,A2,i3,A2,d18.9)
31 format(6x,A6,i3,A2,i3)
write(id_out,'(2X,A)') "Module dip_param"
write(id_out,'(5X,A)') "IMPLICIT NONE"
write(id_out,'(5X,A,I0)') "Integer,parameter :: np=",npar
write(id_out,'(5X,A,I0,A)') "Double precision :: p(",npar,")"
write(id_out,'(5X,A,I0,A)') "integer :: pst(2,",maxpar_keys,")"
write(id_out,'(5X,A)') "contains"
write(id_out,*)''
write (id_out,'(5x,a)') "SUBROUTINE init_dip_planar_data()"
write (id_out,'(8X,A)') "implicit none"
do i=1,npar
write(id_out,30) 'p(',i,')=',p(i)
enddo
do i=1,maxpar_keys
write(id_out,31) 'pst(1,',i,')=',pst(1,i)
write(id_out,31) 'pst(2,',i,')=',pst(2,i)
enddo
write(id_out,"(A)") "End SUBROUTINE init_dip_planar_data"
write(id_out,"(A)") "End Module dip_param"
close(id_out)
end subroutine
!----------------------------------------------------------------------------------------------------
! <subroutine print_ErrorSummary: calculates the rms errros and prints them in the corresponding file
subroutine print_ErrorSummary(id_out,y,ymod,wt)
use dim_parameter,only: nstat,rms_thr,ntot,numdatpt
use io_parameters,only: llen
implicit none
! IN: variables
integer id_out
double precision y(ntot,numdatpt),ymod(ntot,numdatpt)
double precision wt(ntot,numdatpt)
! INTERNAL: variables
! Counter and RMS variables
double precision Cut_thr(nstat)
double precision Output_rms(ntot),Cut_rms(nstat),Weighted_rms
integer Output_num(ntot),Cut_num(nstat)
double precision Weighted_wt
double precision Total_rms,Total_Weighted_rms
double precision Total_Energie_rms,Total_State_rms(nstat)
double precision Cut_Energie_rms, Cut_State_rms(nstat)
! Variables for computing the NRMSE
!double precision:: ymean(ntot),ysum(ntot),NRMSE
! loop control
integer j,pt
! Fabian
character(len=llen) fmt
! initialize RMS variables
Output_rms(1:ntot) = 0.d0
Output_num(1:ntot) = 0
Weighted_rms = 0.d0
Weighted_wt = 0.d0
Cut_rms(1:nstat)= 0.d0
Cut_num(1:nstat)= 0
! Define Threshold for Cut_* RMS Values
Cut_thr(1:nstat) = rms_thr(1:nstat)
! SUMM!
! Loop over all Datapoints
do pt=1,numdatpt
! get unweighted rms for each output value and count their number
do j=1,ntot
if(wt(j,pt).gt.0.d0) then
Output_rms(j) = Output_rms(j) +
> (ymod(j,pt)-y(j,pt))**2
Output_num(j)=Output_num(j) + 1
endif
enddo
! get the unweighted rms under the given threshold and count their number
do j=1,nstat
if(wt(j,pt).gt.0.d0) then
if(y(j,pt).le.Cut_thr(j)) then
Cut_rms(j) = Cut_rms(j) +
> (ymod(j,pt)-y(j,pt))**2
Cut_num(j) = Cut_num(j) + 1
endif
endif
enddo
! get the weighted rms over all output values
Weighted_rms = Weighted_rms +
> sum(((ymod(1:ntot,pt)-y(1:ntot,pt))**2)
> *(wt(1:ntot,pt)**2))
Weighted_wt = Weighted_wt + sum(wt(1:ntot,pt)**2)
enddo
! NORM!
! TOTAL RMS:
! unweighted
Total_rms =
> dsqrt(sum(Output_rms(1:ntot)) /(sum(Output_num(1:ntot))))
! Weighted
Total_Weighted_rms = dsqrt(Weighted_rms/Weighted_wt)
! unweighted, considering only first nstat values
Total_Energie_rms =
> dsqrt(sum(Output_rms(1:nstat)) /(sum(Output_num(1:nstat))))
! unweighted,for each of the first nstat values separatly
Total_State_rms(1:nstat) =
> dsqrt(Output_rms(1:nstat) / Output_num(1:nstat))
! unweighted,first nstat values only counting points under given threshold
Cut_Energie_rms =
> dsqrt(sum(Cut_rms(1:nstat)) /(sum(Cut_num(1:nstat))))
! unweighted,each nstat values seperatly only counting points under threshold
Cut_State_rms(1:nstat) =
> dsqrt(Cut_rms(1:nstat)/Cut_num(1:nstat))
! WRITE!
! make the actual writing into the file
write(id_out,39)
write(id_out,40)
write(id_out,41) Total_rms, Total_rms*au2Debye!Total_rms*h2icm
write(id_out,42) sum(Output_num(1:ntot))
write(id_out,43) Total_Weighted_rms, Total_Weighted_rms*h2icm
write(id_out,44) Weighted_wt
write(id_out,45) Total_Energie_rms, Total_Energie_rms*h2icm
write(id_out,42) sum(Output_num(1:nstat))
write(fmt,'("(A,10x,A,",I2,"f8.1)")') nstat
write(id_out,fmt) '#','State resolved RMS(icm): ',
$ Total_State_rms(1:nstat)*h2icm
write(fmt,'("(A,10x,A,",I2,"i8)")') nstat
write(id_out,fmt) '#','No. of Points per State: ',
$ Output_num(1:nstat)
write(id_out,51)
! write the errors under a given threshold if there were any points
if(any(Cut_num(1:nstat).gt.0)) then
write(id_out,48) Cut_Energie_rms, Cut_Energie_rms*h2icm
write(id_out,42) sum(Cut_num(1:nstat))
write(fmt,'("(A,10x,A,",I2,"f8.1,A)")') nstat
write(id_out,fmt) '#','Red. State resolved RMS: ',
$ Cut_State_rms(1:nstat)*h2icm,' icm'
write(fmt,'("(A,10x,A,",I2,"i8)")') nstat
write(id_out,fmt) '#','No. of Points per State: ',
$ Cut_num(1:nstat)
write(fmt,'("(A,10x,A,",I2,"f8.1,A)")') nstat
write(id_out,fmt) '#','Threshold per State: ',
$ Cut_thr(1:nstat)*h2icm,' icm above Reference Point.'
endif
write(id_out,39)
! FORMAT! specifications for the writing
39 format(250('#'))
40 format('#',10x,'ERROR SUMMARY: ')
41 format('#',10x,'Total RMS: ',g16.8, '(',g16.8,
> ' Debye)')
42 format('#',10x,'No. of Points: ',i10)
43 format('#',10x,'Total weighted RMS: ',g16.8, '(',f8.1,' icm)')
44 format('#',10x,'Sum of point weights: ',f16.8)
45 format('#',10x,'Total Energie RMS: ',g16.8, '(',f8.1,' icm)')
48 format('#',10x,'Red. Energie RMS: ',g16.8,'(',f8.1,' icm)')
51 format('#')
end subroutine
!----------------------------------------------------------------------------------------------------
subroutine print_plotfiles(x,y,wt,ymod)
use dim_parameter,only: ndata,sets,qn,ntot,numdatpt,plot_coord
implicit none
! IN: variables
double precision x(qn,numdatpt),y(ntot,numdatpt)
double precision wt(ntot,numdatpt), ymod(ntot,numdatpt)
! INTERNAL: variables
integer sstart,ssend,set,id_plot
! Initialize position pointer
ssend=0
! loop over datasets and print the plotfiles
do set=1 ,sets
if(ndata(set).eq.0) cycle
id_plot=50+set
call open_plotfile(id_plot,set)
write(id_plot,'(A)') '# -*- truncate-lines: t -*-'
! get start and end point of each set
sstart=ssend+1
ssend=ssend+ndata(set)
if (plot_coord(set).eq.0) then
call print_plotwalk(x(:,sstart:ssend),y(:,sstart:ssend),
> wt(:,sstart:ssend),ymod(:,sstart:ssend),
> ndata(set),id_plot,set)
else
call print_plotcoord(plot_coord(set),
> x(:,sstart:ssend),y(:,sstart:ssend),
> wt(:,sstart:ssend),ymod(:,sstart:ssend),
> ndata(set),id_plot,set)
endif
close(id_plot)
enddo
end subroutine
!----------------------------------------------------------------------------------------------------
subroutine print_plotwalk(x,y,wt,ymod,npt,id_plot,set)
use dim_parameter,only: qn,ntot
use io_parameters,only: llen
implicit none
! IN: variables
integer id_plot,npt,set
double precision x(qn,npt),y(ntot,npt),ymod(ntot,npt),wt(ntot,npt)
! INTERNAL: variables
double precision xdiff(qn),walktime
double precision walknorm
! loop control
integer i,j
character(len=llen) fmt
j=ntot-1
call print_plotheader(id_plot,0,npt,set)
call getwalknorm(x,walknorm,npt)
walktime = 0.d0
do i=1,npt
if(i.gt.1) then
xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1)
walktime = walktime + dsqrt(sum(xdiff(1:qn)**2))/walknorm
endif
write(id_plot,"(ES16.8,*(3(ES16.8),:))")
> walktime ,ymod(:,i),y(:,i),(wt(:,i))
enddo
end subroutine
!----------------------------------------------------------------------------------------------------
subroutine print_plotcoord(coord,x,y,wt,ymod,npt,id_plot,set)
use dim_parameter,only: qn,ntot
use io_parameters,only: llen
implicit none
! IN: variables
integer, intent(in) :: id_plot,npt,set,coord
double precision, intent(in) :: x(qn,npt),y(ntot,npt)
double precision, intent(in) :: ymod(ntot,npt),wt(ntot,npt)
! loop control
integer i
call print_plotheader(id_plot,coord,npt,set)
do i=1,npt
! write(id_plot,"(ES16.8,*(3(ES16.8),:))")
! > x(coord,i), ymod(:,i),y(:,i),(wt(:,i))
write(id_plot,"(2ES16.8,*(3(ES16.8),:))")
> x(coord,i), x(coord+1,i),y(:,i)
enddo
end subroutine
!----------------------------------------------------------------------------------------------------
subroutine print_plotheader(id_plot,coord,npt,set)
use dim_parameter,only: qn,ntot
use io_parameters,only: llen
implicit none
integer, intent(in) :: id_plot,npt,set,coord
character(len=llen) fmt
write(id_plot,'("#SET: ",i5)') set
write(id_plot,'("#OUTPUT VALUES",i4)') ntot
write(id_plot,'("#DATA POINTS: ",i4)') npt
if (coord.le.0) then
write(id_plot,'("#t(x) = WALK")')
else
write(id_plot,'("#t(x) = x(",I0,")")') coord
endif
write(id_plot,'("#UNIT: hartree")')
write(id_plot,'()')
write(id_plot,'("#",A15)',advance='no') "t(x)"
write(fmt,'("(3(7X,A9,",I3,"(16x)))")') ntot-1
write(id_plot,fmt) 'ymod(p,x)','y(x) ','wt(x) '
end subroutine
!----------------------------------------------------------------------------------------------------
! <subroutine walknorm calulates the distance in coordinate space for each set
subroutine getwalknorm(x,walknorm,npt)
use dim_parameter,only: qn
implicit none
! IN: variables
integer npt
double precision x(qn,npt)
double precision walknorm
! INTERNAL: variables
double precision xdiff(qn)
integer i
walknorm =0.d0
do i=2,npt
xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1)
walknorm = walknorm + dsqrt(sum(xdiff(1:qn)**2))
enddo
end subroutine
!----------------------------------------------------------------------------------------------------
! <Subroutine for generating output filenames and openeing the correspondign files
subroutine open_plotfile(id_plot,set)
implicit none
! IN: Variables
integer id_plot,set
! INTERNAL: Variables
character(len=30) name !name of output file
! define name sheme for plot files
if (set .lt. 10 ) then
write(name,203) set
else
write(name,202) set
endif
202 format('scan',I2,'.dat')
203 format('scan0',I1,'.dat')
!write (name,202) set
c open plotfile
open(id_plot,file=name)
end subroutine
!----------------------------------------------------------------------------------------------------
! <Subroutine for generating output filenames and openeing the correspondign files
subroutine open_outfile(id_out,it_makro)
implicit none
integer id_out,it_makro
character(len=30) outname !name of output file
543 format('mnlfit-',i1,'.out')
544 format('mnlfit-',i2,'.out')
545 format('mnlfit-',i3,'.out')
if(it_makro.lt.10) then
write(outname,543) it_makro
else if (it_makro.lt.100) then
write(outname,544) it_makro
else if (it_makro.lt.1000) then
write(outname,545) it_makro
else
write(6,*)
> 'ERROR: No rule for Outputfile naming for MAXIT >= 1000'
stop
endif
open (id_out,file=outname)
end subroutine
!----------------------------------------------------------------------------------------------------
! <Subroutine for printing the Parameterkeys for use in Input File
! < prints the keystring given in keys.incl and the corresponding parameters when there was atleast one parameter given in the input for the spcific key
! < how many parameters and spreads per line are printed can be specified with the hardcoded parameters np and nsp but they must be atleast >=2
! <@param id_out specifies the file in which the Parameters are Printed
! <@param p vector containing one set of parameter values
! <@param p_act vector containing the active state 0 (inactive) or 1 (active) for each parameter
! <@param p_spread vector containing the spreads for each parameter
! <@param npar lenght of the parmeter vectors (p,p_act,p_spread)
! <@TODO extract subroutine for printing the multiline values, would make this more readable
subroutine print_parameterblock(id_out,p,p_act,p_spread,npar)
use dim_parameter,only: pst, facspread
use io_parameters,only: key, parkeynum,parkeylen,llen
implicit none
! IN: Variables
integer id_out,npar,p_act(npar)
double precision p(npar),p_spread(npar)
! INTERNAL: variables
! loop index
integer i,k,l,t,n !< internal variables for loops and positions in parameter vectors
! number of values per line, values must be atleast 2 set this to personal preference
integer, parameter :: np=5,nsp=5
character(len=llen) fmt
! Write header for Parameter block
1 format('!',200('='))
write(id_out,1)
write(id_out,'(A2,5x,A11,i3)') '! ','PARAMETER: ',npar
write(id_out,1)
! loop over all Parameter Keys
do i = 1, parkeynum
! save start and end of parameter block for specific key
k = pst(1,i)
l = pst(1,i)+pst(2,i)-1
! print only used keys with atleast one parameter
if(pst(2,i).gt.0) then
write(fmt,'("(a",I3,"'' ''i3)")') parkeylen
write(id_out,fmt) adjustl(key(1,i)), pst(2,i)
! write the actual parameters -> subroutine print_parameterlines()?
if(l-k.le.(np-1)) then
write(fmt,'("(a",I3,"'' ''",I3,"g24.15)")') parkeylen,np
write(id_out,fmt) key(2,i),(p(n), n=k,l)
else
! start of multi line parameter print, number of values per line specified by np
write(fmt,'("(a",I3,"'' ''",I3,"g24.15'' &'')")')
$ parkeylen,np
write(id_out,fmt) key(2,i),(p(n), n=k,k+(np-1))
t=k+np
! write continuation lines till left parameters fit on last line
do while(t.le.l)
if(l-t.le.(np-1)) then
write(fmt,'("(",I3,"x'' ''",I3,"g24.15)")')
$ parkeylen,np
write(id_out,fmt) (p(n), n=t, l)
else
write(fmt,'("(",I3,"x'' ''",I3,"g24.15'' &'')")')
$ parkeylen,np
write(id_out,fmt) (p(n), n=t, t+(np-1))
endif
t=t+np
enddo
endif !-> end subroutine print_parameterlines
! write parameter active state in one line
write(fmt,'("(a",I3,"'' ''","50i3)")') parkeylen
write(id_out,fmt) key(3,i),(p_act(n),n=k,l)
! write the spreads for each parameter
if(l-k.le.(np-1)) then
write(fmt,'("(a",I3,"'' ''",I3,"g24.8)")') parkeylen,nsp
write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,l)
else
! start of multiline spread values
write(fmt,'("(a",I3,"'' ''",I3,"g24.8'' &'')")')
$ parkeylen,nsp
write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,k
> +(np-1))
t=k+nsp
! write continuation lines till left spreads fit on last line
do while(t.le.l)
if(l-t.le.(np-1)) then
write(fmt,'("(",I3,"x'' ''",I3,"g24.8)")')
$ parkeylen,nsp
write(id_out,fmt) (p_spread(n)/facspread, n=t, l)
else
write(fmt,'("(",I3,"x'' ''",I3,"g24.8'' &'')")')
$ parkeylen,nsp
write(id_out,fmt) (p_spread(n)/facspread, n=t, t
> +(np-1))
endif
t=t+np
enddo
endif
! print empty line between diffrent parameter blocks for better readability
write(id_out,'(" ")')
endif
enddo
end subroutine
!----------------------------------------------------------------------------------------------------
! <Subroutine for printing the current Parameters and their active state
! < prints only the numeric values of the parameters and does not specify the corresponding key
! <@param npar number of parameter
! <@param id_out specifies the output file
! <@param p,p_act parameter vectors containing the values and the activity state of parameters
subroutine print_parameterstate(id_out,p,p_act,npar)
implicit none
! IN: Variables
integer npar,id_out
double precision p(npar)
integer p_act(npar)
! INTERNAL: Variables
integer i !< loop control
integer nopt !< number of counted active parameters
character(len=16) opt(npar) !< string for optimisation state
! initialize number of opt parameters and the string vector opt
nopt=0
opt = ' not opt. '
! loop over all parameters and check their active state count if active and set string to opt
do i=1,npar
! Nicole: change due to value 2 of p_act
! if(p_act(i).eq.1) then
if(p_act(i).ge.1) then
opt(i) = ' opt. '
nopt=nopt+1
endif
enddo
! print the Parameters and their active state within separating lines
write(id_out,*)''
write(id_out,block_line)
write(id_out,*) 'Parameters:'
write(id_out,sep_line)
write(id_out,'(5g14.6)') (p(i),i=1,npar)
write(id_out,'(5a14)') (opt(i),i=1,npar)
write(id_out,sep_line)
write(id_out,'("No. of optimized parameters: ",i6)') nopt
write(id_out,block_line)
write(id_out,*)''
end subroutine
!----------------------------------------------------------------------------------------------------
! <Subroutine for printing coordinates,refdata,modeldata,diffrence between them and the weights
! <@param id_out identiefies the output file
! <@param x vector of input pattern for each datapoint
! <@param y vector of expected output patterns for each datapoint
! <@param ymod vector of output patterns generated by the model depending on paramerters
! <@param wt vector of weights for each datapoint
! <@param qn number of input patterns
! <@param ntot total number of output patterns for each datapoint
! <@param numdatpt number of totatl datapoints
! <@param sets number of sets the datapoints are divided into
! <@param ndata vector containing the number of included datapoints for each set
! <@param i,j,point internal variables for loop controll and datapoint counting
subroutine print_data(id_out,x,y,ymod,wt)
use dim_parameter,only: sets,ndata,qn,ntot,numdatpt,qn_read
implicit none
! IN: Variables
integer id_out
double precision x(qn,numdatpt)
double precision y(ntot,numdatpt),ymod(ntot,numdatpt)
double precision wt(ntot,numdatpt)
! INTERNAL: Variables
integer i,j,point
18 format(A8,i6)
19 format (3(A15,3x), 2x, A18 , 4x, A12)
! print seperating line and header for Data output
write(id_out,*) 'Printing Data Sets:'
write(id_out,19) adjustl('y(x)'),adjustl('ymod(x)'),
> adjustl('y(x)-ymod(x)'),adjustl('weight'),
> adjustl('x(1:qn_read) ')
write(id_out,sep_line)
! loop over all datapoints for each set and count the actual datapointnumber with point
point=0
do i=1,sets
write(id_out,18) 'Set: ', i
do j=1,ndata(i)
write(id_out,18) 'Point: ', j
point=point+1
! print all data for one datapoint
call print_datapoint(id_out,x(:,point),y(:,point),
> ymod(:,point),wt(:,point))
write(id_out,sep_line)
enddo
enddo
! write end of data statement and two seperating lines
write(id_out,block_line)
write(id_out,*) ''
end subroutine
!----------------------------------------------------------------------------------------------------
! <Subroutine prints a single Datapoint splits Data in nstat nci(ndiab) blocks for readability
! <@param id_out identiefies the output file
! <@param x vector of input pattern for each datapoint
! <@param y vector of expected output patterns for each datapoint
! <@param ymod vector of output patterns generated by the model depending on paramerters
! <@param wt vector of weights for each datapoint
! <@param qn number of input patterns
! <@param ntot total number of output patterns for each datapoint
! <@param i,j,k internal variables for loop controll and counting
subroutine print_datapoint(id_out,x,y,ymod,wt)
use dim_parameter,only: nstat,ndiab,nci,qn,ntot,qn_read
use io_parameters,only: llen
implicit none
integer id_out
double precision x(qn),y(ntot),ymod(ntot),wt(ntot)
integer i,j,k
18 format(A10,i3)
19 format(3F18.8, 2X, F18.6, 4X,*(F12.6))
! print the nstat output patterns
do i=1,nstat
write(id_out,19)y(i),ymod(i),ymod(i)-y(i), wt(i), x(1:qn)
enddo
! loop over number (nci) of metadata with lenght (ndiab)
do i=1,nci
write(id_out,18) 'nci: ',i
do j=1,ndiab
k=nstat + (i-1)*ndiab + j
write(id_out,19) y(k),ymod(k),(ymod(k)-y(k)),
> wt(k), x(1:qn_read)
enddo
enddo
end subroutine
end module write_mod