Compare commits
3 Commits
main
...
model_muxy
| Author | SHA1 | Date |
|---|---|---|
jean paul nshuti
|
412bb832f2 | |
|
jean paul nshuti
|
464dca1bbb | |
|
jean paul nshuti
|
733a5d4172 |
|
|
@ -0,0 +1,164 @@
|
|||
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
|
||||
#VPATH += $(src)model/PES
|
||||
######################################################################
|
||||
|
||||
#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 matrix_form.o ctrans.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) )
|
||||
|
||||
|
||||
#Note: Since we are using modules, you have carefully choose the order of compilation and take dependencies between modules and subroutines into account!
|
||||
|
||||
######################################################################
|
||||
|
||||
# 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)
|
||||
$(FC) $(FFLAGS) $(random_obj) $(objects) -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)
|
||||
|
||||
all: clean $(main) umb_genetic
|
||||
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,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=[1.0228710942d0,120.d0])
|
||||
|
||||
!##########################################################################
|
||||
! coordinate order; the first #I number of coords are given to the
|
||||
! ANN, where #I is the number of input neurons. The position i in
|
||||
! pat_index corresponds to a coordinate, the value of pat_index(i)
|
||||
! signifies its position.
|
||||
!
|
||||
! The vector is ordered as follows:
|
||||
! a,xs,ys,xb,yb,b,rs**2,rb**2,b**2,
|
||||
! es*eb, es**3, eb**3,es**2*eb, es*eb**2
|
||||
! ri**2 := xi**2+yi**2 = ei**2; ei := (xi,yi), i = s,b
|
||||
!
|
||||
! parts not supposed to be read by ANN are marked by ';' for your
|
||||
! convenience.
|
||||
!##########################################################################
|
||||
! a,rs**2,rb**2,es*eb,es**3,eb**3,es**2*eb,es*eb**2,b**2 #I=9 (6D)
|
||||
parameter (pat_index=[1,2,3,4,5,6,7,8,9,10,11,12,13,14])
|
||||
!**************************************************************************
|
||||
|
|
@ -0,0 +1,98 @@
|
|||
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 diabmodel, 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.)
|
||||
! 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
|
||||
double precision:: alpha
|
||||
integer:: lda,ldb,beta,ldc
|
||||
double precision,dimension(ndiab,ndiab):: temp1,temp2
|
||||
call diab(ex,ey,n,x1_m(:,n),x2_m(:,n),p)
|
||||
! init eigenvector matrix
|
||||
|
||||
u = 0.d0
|
||||
skip=.false.
|
||||
ymod=0.0d0
|
||||
call Full_diab_upper(ex,ey,ymod)
|
||||
!write(*,'(16f10.4)') ex
|
||||
! ymod(1)=ex(1,1)
|
||||
! ymod(2)=ex(1,2)
|
||||
! ymod(3)=ex(1,3)
|
||||
! ymod(4)=ex(1,4)
|
||||
! ymod(5)=ex(2,2)
|
||||
! ymod(6)=ex(2,3)
|
||||
! ymod(7)=ex(2,4)
|
||||
! ymod(8)=ex(3,3)
|
||||
! ymod(9)=ex(3,4)
|
||||
! ymod(10)=ex(4,4)
|
||||
|
||||
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
|
||||
|
|
@ -0,0 +1,367 @@
|
|||
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)
|
||||
cart(:,1)=0.0d0
|
||||
cart(1:3,2:4) = reshape([ q(1:9) ], shape(cart(1:3,2:4)))
|
||||
call cart2int(cart,qint)
|
||||
x1(1:qn)=qint(1:qn)
|
||||
!x1(6)=-x1(6)
|
||||
x2(1:qn)=0.0d0 !qint(1:qn)
|
||||
!x1=q
|
||||
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
|
||||
|
||||
|
|
@ -0,0 +1,40 @@
|
|||
! <subroutine for manipulating the input Data before the Fit
|
||||
subroutine data_transform(q,x1,x2,y,wt,p,npar,p_act)
|
||||
use dim_parameter,only : nstat,pst,ntot,qn,numdatpt,ndiab
|
||||
use ctrans_mod, only: ctrans
|
||||
use data_matrix
|
||||
! use david_ctrans_mod, only: ctrans_d
|
||||
implicit none
|
||||
! IN: variables
|
||||
integer npar
|
||||
double precision q(qn,numdatpt),x1(qn,numdatpt),x2(qn,numdatpt)
|
||||
double precision y(ntot,numdatpt),wt(ntot,numdatpt)
|
||||
double precision p(npar),mat_x(ndiab,ndiab),mat_y(ndiab,ndiab)
|
||||
double precision v(ndiab,ndiab),E(nstat)
|
||||
integer p_act(npar), pt
|
||||
logical dbg
|
||||
parameter (dbg=.false.)
|
||||
integer TYPES,BLK ! TYPE OF THE CALCULATION AND THE BLOCK IF TYEPE IS 3
|
||||
double precision U(ndiab,ndiab), U_ref(ndiab,ndiab) ! Transformation matrix
|
||||
|
||||
|
||||
|
||||
! get the ref transformation matrix
|
||||
!call eval_surface(E,V,U_ref,q(1:qn,1),p)
|
||||
|
||||
do pt=1,numdatpt
|
||||
call ctrans(q(1:qn,pt),x1(:,pt),x2(:,pt))! ctrans the dipole cooordinate.
|
||||
!call ctrans_pes(q(1:qn,pt),x1(:,pt),x2(:,pt))
|
||||
write(7,'(6f18.8)') x1(1:6,pt)
|
||||
y(11:ntot,pt)=-y(11:ntot,pt)
|
||||
enddo
|
||||
|
||||
call weight(wt,y)
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,88 @@
|
|||
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,108)
|
||||
|
||||
character(len=16) parname(24)
|
||||
integer i,j
|
||||
! Defining keys for potential
|
||||
! the electronic states are ordered as: A2" E' and A1'
|
||||
! the name convention here is : A2 E1 A1
|
||||
|
||||
! Naming
|
||||
!--------------------
|
||||
! V: V-TERM OR diagonal term
|
||||
! J: Jahn teller coupling term in E
|
||||
! P: pseudo jahn teller between As and E
|
||||
|
||||
! S: it involves the symmetric term of x**2+y**2
|
||||
! N: It does not involve symmetric term
|
||||
|
||||
|
||||
! diagonal term for 4 states
|
||||
! no zeroth order in V
|
||||
|
||||
|
||||
parname( 1)='VA2N1' ! order 1
|
||||
parname( 2)='VE1N1' ! order 1 witH N
|
||||
parname( 3)='VA1N1' ! order 1 with N
|
||||
parname( 4)='VA2N2' ! order 2
|
||||
parname( 5)='VE1N2' ! order 2 witH N
|
||||
parname( 6)='VA1N2' ! order 2 with N
|
||||
parname( 7)='VA2N3' ! order 3
|
||||
parname( 8)='VE1N3' ! order 3 witH N
|
||||
parname( 9)='VA1N3' ! order 3 with N
|
||||
! Jahn teller within E
|
||||
parname(10)='JE1N0' ! order 0 with N
|
||||
parname(11)='JE1N1' ! order 1 with N
|
||||
parname(12)='JE1N2' ! order 2
|
||||
parname(13)='JE1N3' ! order 3 ! this has 8 terms
|
||||
|
||||
|
||||
! PSeudo Jahn teller couplings
|
||||
|
||||
! coupling of A2 with E
|
||||
parname(14)='PA2E1N0' ! order 0 ! is not there
|
||||
parname(15)='PA2E1N1' ! order 1
|
||||
parname(16)='PA2E1N2' ! Order 2
|
||||
parname(17)='PA2E1N3' !
|
||||
! coupling of A2 with A1
|
||||
!parname(17)='PA2A1N0' ! order 0 ! is not there
|
||||
parname(18)='PA2A1N1' ! order 1
|
||||
parname(19)='PA2A1N2' ! order 2
|
||||
! no order 3 for A2 with A1
|
||||
! coupling of A1 with other
|
||||
! A2 with A1 is already included above
|
||||
|
||||
parname(20)='PE1A1N0' ! order 0
|
||||
parname(21)='PE1A1N1' ! order 1
|
||||
parname(22)='PE1A1N2' ! order 2
|
||||
parname(23)='PE1A1N3' ! order 3
|
||||
|
||||
|
||||
parname(24)='TYPE_CAL'
|
||||
|
||||
do i=1,22
|
||||
do j=1,4
|
||||
key(j, i)=prefix(j)//trim(parname(i))//':' ! first 86 keys are the potential keys
|
||||
!write(*,*) key(j, i)
|
||||
enddo
|
||||
!write(*,*) ''
|
||||
enddo
|
||||
!do i=1,108
|
||||
! do j=1,4
|
||||
! write(*,*) key(j,i)
|
||||
! enddo
|
||||
! write(*,*) ""
|
||||
! enddo
|
||||
end subroutine
|
||||
|
||||
end module keys_mod
|
||||
|
|
@ -0,0 +1,301 @@
|
|||
module data_matrix
|
||||
use dim_parameter, only:ndiab,nstat,ntot,pst,qn
|
||||
! 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
|
||||
|
||||
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)
|
||||
IMPLICIT NONE
|
||||
integer:: ii,i,j
|
||||
double precision, intent(in):: y(:)
|
||||
double precision,intent(out):: Mx(ndiab,ndiab),My(ndiab,ndiab)
|
||||
|
||||
!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=i,ndiab
|
||||
! !mx
|
||||
|
||||
mx(i,j)=y(ii)
|
||||
! ! My
|
||||
my(i,j)=y((ntot/2)+ii)
|
||||
!
|
||||
ii=ii+1
|
||||
enddo
|
||||
enddo
|
||||
call coppy_2_low(mx)
|
||||
call coppy_2_low(my)
|
||||
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 adiabatic_transform(mx,my,U)
|
||||
implicit none
|
||||
double precision, intent(inout) :: mx(ndiab,ndiab), my(ndiab,ndiab)
|
||||
double precision, dimension(:,:), intent(inout) :: U
|
||||
double precision, dimension(ndiab,ndiab) :: temp1, temp2
|
||||
integer i, j
|
||||
! 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(in) :: mx,my
|
||||
double precision,dimension(:),intent(out) :: Yres
|
||||
double precision,dimension(size(mx,1),size(mx,2)) :: 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)
|
||||
Yres = 0.0d0
|
||||
Allocate(mux,source=mx)
|
||||
call DSYEV('V', 'U', size(mx,1), mux, size(mx,1), vx, work, lwork, info)
|
||||
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),1)
|
||||
Yres(size(mx,1)+1:2*size(mx,1)) = vy(1:size(my,1),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,blk)
|
||||
implicit none
|
||||
double precision, intent(inout):: Y(:)
|
||||
double precision, intent(in) :: mx(ndiab,ndiab), my(ndiab,ndiab)
|
||||
integer, intent(in) :: blk
|
||||
integer i,j,ii,nn
|
||||
y=0.0d0
|
||||
select case (blk)
|
||||
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 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(3)
|
||||
! Filling the A1 coupling with E2
|
||||
y(1)=mx(2,4)
|
||||
y(2)=mx(3,4)
|
||||
y(3)=mx(4,2)
|
||||
y(4)=mx(4,3)
|
||||
! my
|
||||
y(5)=my(2,4)
|
||||
y(6)=my(3,4)
|
||||
y(7)=my(4,2)
|
||||
y(8)=my(4,3)
|
||||
case(4)
|
||||
! filling the block of A2 coupling with Es
|
||||
y(1)=mx(1,1)
|
||||
y(2)=mx(1,4)
|
||||
y(3)=mx(4,4)
|
||||
! my
|
||||
y(5)=my(1,1)
|
||||
y(6)=my(1,4)
|
||||
y(7)=my(4,4)
|
||||
|
||||
|
||||
case default
|
||||
write(*,*) "Error in block_diab subroutine, block not recognized"
|
||||
write(*,*) "The block is:", blk
|
||||
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,q)
|
||||
implicit none
|
||||
double precision, intent(inout) :: U(ndiab,ndiab)
|
||||
double precision, intent(in) :: q(qn)
|
||||
integer i,max_row
|
||||
logical,parameter :: dbg_sign =.true.
|
||||
double precision :: theta
|
||||
do i=1,ndiab
|
||||
max_row = maxloc(abs(U(:,i)),1)
|
||||
if (U(max_row,i) .lt. 0) then
|
||||
U(:,i) = -1.0*U(:,i)
|
||||
endif
|
||||
enddo
|
||||
if (dbg_sign) then
|
||||
theta=atan(q(3)/q(2))
|
||||
U=sign(1.0d0,cos(theta))*sign(1.0d0,sin(theta))*U
|
||||
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,108)))
|
||||
blk = int(p(pst(1,108)+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: E' BLOCK"
|
||||
ELSE IF (BLK ==2) THEN
|
||||
write(id_write,*) "Type of calculation: COUPLING BLOCK 1 & 2"
|
||||
ELSE IF (BLK ==3) THEN
|
||||
write(id_write,*) "Type of calculation: cOUPLING BLOCK 4 &2 "
|
||||
ELSE IF (BLK ==4) THEN
|
||||
write(id_write,*) "Type of calculation: A(11),A(14),A(44) ELEMENTS"
|
||||
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)
|
||||
|
||||
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(17:20)=v(1:4) ! copy the first 4 elements of v to y
|
||||
end subroutine
|
||||
|
||||
|
||||
|
||||
end module
|
||||
|
|
@ -0,0 +1,375 @@
|
|||
module diabmodel
|
||||
use dim_parameter,only:qn,ndiab,pst
|
||||
use accuracy_constants, only:dp,idp
|
||||
implicit none
|
||||
logical :: debug=.false.
|
||||
contains
|
||||
|
||||
subroutine diab(ex,ey,n,x1,x2,p)
|
||||
use ctrans_mod, only:ctrans
|
||||
integer,intent(in),optional :: n ! number of parameter & nmbr of points \
|
||||
integer id
|
||||
integer key,i,j
|
||||
double precision, intent(in)::x1(qn),x2(qn)
|
||||
double precision, contiguous,intent(in):: p(:)! array containing parameters
|
||||
double precision, intent(out)::ex(ndiab,ndiab),ey(ndiab,ndiab)
|
||||
key =1
|
||||
call diab_x(ex,x1,x2,key,p)
|
||||
!ey=0.0d0
|
||||
call diab_y(ey,x1,x2,key,p)
|
||||
end subroutine
|
||||
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
subroutine diab_x(e,q,t,key,p)
|
||||
real(dp),intent(in)::q(qn),t(qn)
|
||||
real(dp),intent(out)::e(:,:)
|
||||
integer(idp),intent(in)::key
|
||||
real(dp),intent(in),contiguous::p(:)
|
||||
integer(idp) id,i,j
|
||||
real(dp) tmp_v,xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
|
||||
xs=q(2)
|
||||
ys=q(3)
|
||||
xb=q(4)
|
||||
yb=q(5)
|
||||
a=q(1)
|
||||
b=q(6)
|
||||
|
||||
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=key !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,t,key,p)
|
||||
!integer(idp), intent(in)::npar
|
||||
real(dp),intent(in)::q(qn),t(qn)
|
||||
real(dp),intent(out)::e(:,:)
|
||||
integer(idp),intent(in):: key
|
||||
real(dp),intent(in),contiguous::p(:)
|
||||
integer(idp) id,i,j
|
||||
real(dp) tmp_v,ys,xb,a,b,xs,yb,ss,sb,v3_vec(8)
|
||||
xs=q(2)
|
||||
ys=q(3)
|
||||
xb=q(4)
|
||||
yb=q(5)
|
||||
a=q(1)
|
||||
b=q(6)
|
||||
|
||||
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=key !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
|
||||
|
|
@ -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)
|
||||
|
|
@ -0,0 +1,104 @@
|
|||
module keys_mod
|
||||
implicit none
|
||||
contains
|
||||
subroutine init_keys
|
||||
use io_parameters, only: key
|
||||
character(len=1) prefix(4)
|
||||
parameter (prefix=['N','P','A','S'])
|
||||
!character (len=20) key(4,56)
|
||||
|
||||
character(len=16) parname(56)
|
||||
integer i,j
|
||||
|
||||
! the electronic states are ordered as: A2" E' and A1'
|
||||
! the name convention here is : A2 E1 A1
|
||||
|
||||
! Naming
|
||||
!--------------------
|
||||
! V: V-TERM OR diagonal term
|
||||
! J: Jahn teller coupling term in E
|
||||
! P: pseudo jahn teller between As and E
|
||||
|
||||
! S: it involves the symmetric term of x**2+y**2
|
||||
! N: It does not involve symmetric term
|
||||
|
||||
|
||||
! diagonal term for 4 states
|
||||
!parname( 1)='VA2N0' ! order 0
|
||||
!parname( 2)='VA2S0' ! order 0 with S
|
||||
!parname( 3)='VE1N0' ! order 0 witH N
|
||||
!parname( 4)='VE1S0' ! order 0 with S
|
||||
!parname( 5)='VA1N0' ! order 0 with N
|
||||
!parname( 6)='VA1S0' ! order 0 with S
|
||||
parname( 7)='VA2N1' ! order 1
|
||||
parname( 8)='VA2S1' ! order 1 with S
|
||||
parname( 9)='VE1N1' ! order 1 witH N
|
||||
parname(10)='VE1S1' ! order 1 with S
|
||||
parname(11)='VA1N1' ! order 1 with N
|
||||
parname(12)='VA1S1' ! order 1 with S
|
||||
parname(13)='VA2N2' ! order 2
|
||||
parname(14)='VA2S2' ! order 2 with S ! only 2 term
|
||||
parname(15)='VE1N2' ! order 2 witH N
|
||||
parname(16)='VE1S2' ! order 2 with S
|
||||
parname(17)='VA1N2' ! order 2 with N
|
||||
parname(18)='VA1S2' ! order 2 with S
|
||||
!parname(19)='VA2N3' ! order 3
|
||||
!parname(20)='VA2S3' ! order 3 with S
|
||||
!parname(21)='VE1N3' ! order 3 witH N
|
||||
!parname(22)='VE1S3' ! order 3 with S
|
||||
!parname(23)='VA1N3' ! order 3 with N
|
||||
!parname(24)='VA1S3' ! order 3 with S
|
||||
|
||||
! Jahn teller within E
|
||||
|
||||
parname(25)='JE1N0' ! order 0 with N
|
||||
parname(26)='JE1S0' ! order 0 with S
|
||||
parname(27)='JE1N1' ! order 1 with N
|
||||
parname(28)='JE1S1' ! order 1 with S
|
||||
parname(29)='JE1N2' ! order 2
|
||||
parname(30)='JE1S2' ! order 2
|
||||
parname(31)='JE1N3' ! order 3 ! this has 8 terms
|
||||
!parname(32)='JE1S3' ! order 3 ! i do not have this term
|
||||
|
||||
! PSeudo Jahn teller couplings
|
||||
|
||||
! coupling of A2 with other
|
||||
!parname(33)='PA2E1N0' ! order 0 ! is not there
|
||||
! parname(34)='PA2E1S0' ! order 0 ! is not there
|
||||
! parname(35)='PA2A1N0' ! ORDER 0
|
||||
! parname(36)='PA2A1S0' ! order 0
|
||||
!parname(37)='PA2E1N1' ! order 1
|
||||
!parname(38)='PA2E1S1' ! order 1
|
||||
parname(39)='PA2A1N1' ! order 1
|
||||
parname(40)='PA2A1S1' ! order 1
|
||||
!parname(41)='PA2E1N2' ! order 2
|
||||
!parname(42)='PA2E1S2' ! order 2
|
||||
parname(43)='PA2A1N2'
|
||||
parname(44)='PA2A1S2'
|
||||
parname(45)='PA2E1N3' ! order 3
|
||||
!parname(46)='PA2E1S3' ! order 3
|
||||
parname(47)='PA2A1N3'
|
||||
!parname(48)='PA2A1S3'
|
||||
|
||||
|
||||
|
||||
! coupling of A1 with other
|
||||
! A2 with A1 is already included above
|
||||
|
||||
parname(49)='PE1A1N0' ! order 1
|
||||
parname(50)='PE1A1S0'
|
||||
!parname(51)='PE1A1N1'
|
||||
!parname(52)='PE1A1S1'
|
||||
!parname(53)='PE1A1N2'
|
||||
!parname(54)='PE1A1S2'
|
||||
parname(55)='PE1A1N3'
|
||||
!parname(56)='PE1A1S3'
|
||||
|
||||
do i=1,56
|
||||
do j=1,4
|
||||
key(i, j)=prefix(j)//trim(parname(i))//':'
|
||||
enddo
|
||||
enddo
|
||||
end subroutine
|
||||
|
||||
end module keys_mod
|
||||
|
|
@ -0,0 +1,357 @@
|
|||
module diabmodel
|
||||
use dim_parameter,only:qn,ndiab,pst
|
||||
use accuracy_constants, only:dp,idp
|
||||
implicit none
|
||||
logical :: debug=.false.
|
||||
contains
|
||||
|
||||
subroutine diab(ex,ey,n,x1,x2,p)
|
||||
use ctrans_mod, only:ctrans
|
||||
integer,intent(in),optional :: n ! number of parameter & nmbr of points \
|
||||
integer id
|
||||
integer key,i,j
|
||||
double precision, intent(in)::x1(qn),x2(qn)
|
||||
double precision, contiguous,intent(in):: p(:)! array containing parameters
|
||||
double precision, intent(out)::ex(ndiab,ndiab),ey(ndiab,ndiab)
|
||||
key =87
|
||||
call diab_x(ex,x1,x2,key,p)
|
||||
!ey=0.0d0
|
||||
call diab_y(ey,x1,x2,key,p)
|
||||
end subroutine
|
||||
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
subroutine diab_x(e,q,t,key,p)
|
||||
|
||||
|
||||
real(dp),intent(in)::q(qn),t(qn)
|
||||
real(dp),intent(out)::e(:,:)
|
||||
integer(idp),intent(in)::key
|
||||
real(dp),intent(in),contiguous::p(:)
|
||||
integer(idp) id,i,j
|
||||
real(dp) tmp_v,xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
|
||||
xs=q(2)
|
||||
ys=q(3)
|
||||
xb=q(4)
|
||||
yb=q(5)
|
||||
a=q(1)
|
||||
b=q(6)
|
||||
|
||||
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=key !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 + b**2* &
|
||||
(p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
|
||||
id=id+1 !8
|
||||
e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb+ b**2* &
|
||||
(p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
|
||||
e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb+ b**2* &
|
||||
(p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
|
||||
id=id+1 !9
|
||||
e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + b**2* &
|
||||
(p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
|
||||
|
||||
! 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 + &
|
||||
b**2*(p(pst(1,id)+5) +b**4*(p(pst(1,id)+6)) + b**6*(p(pst(1,id)+7)) + &
|
||||
b**8*(p(pst(1,id)+8)))
|
||||
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) -&
|
||||
b**8*(p(pst(1,id)+5))
|
||||
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))
|
||||
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))
|
||||
|
||||
!! THE COUPLING OF A2 WITH A1
|
||||
!####################################################
|
||||
!####################################################
|
||||
! order 1
|
||||
id=id+1 !17
|
||||
e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
|
||||
id=id+1 !18
|
||||
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 !19
|
||||
e(2,4)=e(2,4)+p(pst(1,id))
|
||||
|
||||
! order 1
|
||||
id=id+1 !20
|
||||
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 !21
|
||||
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)
|
||||
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)
|
||||
|
||||
!! 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,t,key,p)
|
||||
!integer(idp), intent(in)::npar
|
||||
real(dp),intent(in)::q(qn),t(qn)
|
||||
real(dp),intent(out)::e(:,:)
|
||||
integer(idp),intent(in):: key
|
||||
real(dp),intent(in),contiguous::p(:)
|
||||
integer(idp) id,i,j
|
||||
real(dp) tmp_v,ys,xb,a,b,xs,yb,ss,sb,v3_vec(8)
|
||||
xs=q(2)
|
||||
ys=q(3)
|
||||
xb=q(4)
|
||||
yb=q(5)
|
||||
a=q(1)
|
||||
b=q(6)
|
||||
|
||||
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=key !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 +b**2* &
|
||||
(p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
|
||||
id=id+1 !8
|
||||
e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb+b**2* &
|
||||
(p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
|
||||
e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* &
|
||||
(p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
|
||||
id=id+1 !9
|
||||
e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* &
|
||||
(p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
|
||||
|
||||
! 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 +&
|
||||
b**8*(p(pst(1,id)+5))
|
||||
! 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))
|
||||
! 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)
|
||||
! 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
|
||||
|
|
@ -0,0 +1,110 @@
|
|||
NVA2N1:
|
||||
PVA2N1:
|
||||
AVA2N1:
|
||||
SVA2N1:
|
||||
|
||||
NVE1N1:
|
||||
PVE1N1:
|
||||
AVE1N1:
|
||||
SVE1N1:
|
||||
|
||||
NVA1N1:
|
||||
PVA1N1:
|
||||
AVA1N1:
|
||||
SVA1N1:
|
||||
|
||||
NVA2N2:
|
||||
PVA2N2:
|
||||
AVA2N2:
|
||||
SVA2N2:
|
||||
|
||||
NVE1N2:
|
||||
PVE1N2:
|
||||
AVE1N2:
|
||||
SVE1N2:
|
||||
|
||||
NVA1N2:
|
||||
PVA1N2:
|
||||
AVA1N2:
|
||||
SVA1N2:
|
||||
|
||||
NVA2N3:
|
||||
PVA2N3:
|
||||
AVA2N3:
|
||||
SVA2N3:
|
||||
|
||||
NVE1N3:
|
||||
PVE1N3:
|
||||
AVE1N3:
|
||||
SVE1N3:
|
||||
|
||||
NVA1N3:
|
||||
PVA1N3:
|
||||
AVA1N3:
|
||||
SVA1N3:
|
||||
|
||||
NJE1N0:
|
||||
PJE1N0:
|
||||
AJE1N0:
|
||||
SJE1N0:
|
||||
|
||||
NJE1N1:
|
||||
PJE1N1:
|
||||
AJE1N1:
|
||||
SJE1N1:
|
||||
|
||||
NJE1N2:
|
||||
PJE1N2:
|
||||
AJE1N2:
|
||||
SJE1N2:
|
||||
|
||||
NJE1N3:
|
||||
PJE1N3:
|
||||
AJE1N3:
|
||||
SJE1N3:
|
||||
|
||||
NPA2E1N0:
|
||||
PPA2E1N0:
|
||||
APA2E1N0:
|
||||
SPA2E1N0:
|
||||
|
||||
NPA2E1N1:
|
||||
PPA2E1N1:
|
||||
APA2E1N1:
|
||||
SPA2E1N1:
|
||||
|
||||
NPA2E1N2:
|
||||
PPA2E1N2:
|
||||
APA2E1N2:
|
||||
SPA2E1N2:
|
||||
|
||||
NPA2A1N1:
|
||||
PPA2A1N1:
|
||||
APA2A1N1:
|
||||
SPA2A1N1:
|
||||
|
||||
NPA2A1N2:
|
||||
PPA2A1N2:
|
||||
APA2A1N2:
|
||||
SPA2A1N2:
|
||||
|
||||
NPE1A1N0:
|
||||
PPE1A1N0:
|
||||
APE1A1N0:
|
||||
SPE1A1N0:
|
||||
|
||||
NPE1A1N1:
|
||||
PPE1A1N1:
|
||||
APE1A1N1:
|
||||
SPE1A1N1:
|
||||
|
||||
NPE1A1N2:
|
||||
PPE1A1N2:
|
||||
APE1A1N2:
|
||||
SPE1A1N2:
|
||||
|
||||
NTYPE_CAL:
|
||||
PTYPE_CAL:
|
||||
ATYPE_CAL:
|
||||
STYPE_CAL:
|
||||
|
||||
|
|
@ -0,0 +1,540 @@
|
|||
NEXITEN:
|
||||
PEXITEN:
|
||||
AEXITEN:
|
||||
SEXITEN:
|
||||
|
||||
NTMC_CH:
|
||||
PTMC_CH:
|
||||
ATMC_CH:
|
||||
STMC_CH:
|
||||
|
||||
NEVA1:
|
||||
PEVA1:
|
||||
AEVA1:
|
||||
SEVA1:
|
||||
|
||||
NEVU:
|
||||
PEVU:
|
||||
AEVU:
|
||||
SEVU:
|
||||
|
||||
NEVE1:
|
||||
PEVE1:
|
||||
AEVE1:
|
||||
SEVE1:
|
||||
|
||||
NEVE2:
|
||||
PEVE2:
|
||||
AEVE2:
|
||||
SEVE2:
|
||||
|
||||
NEVA1U:
|
||||
PEVA1U:
|
||||
AEVA1U:
|
||||
SEVA1U:
|
||||
|
||||
NEVA1E1:
|
||||
PEVA1E1:
|
||||
AEVA1E1:
|
||||
SEVA1E1:
|
||||
|
||||
NEVA1E2:
|
||||
PEVA1E2:
|
||||
AEVA1E2:
|
||||
SEVA1E2:
|
||||
|
||||
NEVUE1:
|
||||
PEVUE1:
|
||||
AEVUE1:
|
||||
SEVUE1:
|
||||
|
||||
NEVUE2:
|
||||
PEVUE2:
|
||||
AEVUE2:
|
||||
SEVUE2:
|
||||
|
||||
NEVE1E2:
|
||||
PEVE1E2:
|
||||
AEVE1E2:
|
||||
SEVE1E2:
|
||||
|
||||
NEVA1UE1:
|
||||
PEVA1UE1:
|
||||
AEVA1UE1:
|
||||
SEVA1UE1:
|
||||
|
||||
NEVA1UE2:
|
||||
PEVA1UE2:
|
||||
AEVA1UE2:
|
||||
SEVA1UE2:
|
||||
|
||||
NEVA1E1E2:
|
||||
PEVA1E1E2:
|
||||
AEVA1E1E2:
|
||||
SEVA1E1E2:
|
||||
|
||||
NEVUE1E2:
|
||||
PEVUE1E2:
|
||||
AEVUE1E2:
|
||||
SEVUE1E2:
|
||||
|
||||
NEVA1UE1E2:
|
||||
PEVA1UE1E2:
|
||||
AEVA1UE1E2:
|
||||
SEVA1UE1E2:
|
||||
|
||||
NA2VA1:
|
||||
PA2VA1:
|
||||
AA2VA1:
|
||||
SA2VA1:
|
||||
|
||||
NA2VU:
|
||||
PA2VU:
|
||||
AA2VU:
|
||||
SA2VU:
|
||||
|
||||
NA2VE1:
|
||||
PA2VE1:
|
||||
AA2VE1:
|
||||
SA2VE1:
|
||||
|
||||
NA2VE2:
|
||||
PA2VE2:
|
||||
AA2VE2:
|
||||
SA2VE2:
|
||||
|
||||
NA2VA1U:
|
||||
PA2VA1U:
|
||||
AA2VA1U:
|
||||
SA2VA1U:
|
||||
|
||||
NA2VA1E1:
|
||||
PA2VA1E1:
|
||||
AA2VA1E1:
|
||||
SA2VA1E1:
|
||||
|
||||
NA2VA1E2:
|
||||
PA2VA1E2:
|
||||
AA2VA1E2:
|
||||
SA2VA1E2:
|
||||
|
||||
NA2VUE1:
|
||||
PA2VUE1:
|
||||
AA2VUE1:
|
||||
SA2VUE1:
|
||||
|
||||
NA2VUE2:
|
||||
PA2VUE2:
|
||||
AA2VUE2:
|
||||
SA2VUE2:
|
||||
|
||||
NA2VE1E2:
|
||||
PA2VE1E2:
|
||||
AA2VE1E2:
|
||||
SA2VE1E2:
|
||||
|
||||
NA2VA1UE1:
|
||||
PA2VA1UE1:
|
||||
AA2VA1UE1:
|
||||
SA2VA1UE1:
|
||||
|
||||
NA2VA1UE2:
|
||||
PA2VA1UE2:
|
||||
AA2VA1UE2:
|
||||
SA2VA1UE2:
|
||||
|
||||
NA2VA1E1E2:
|
||||
PA2VA1E1E2:
|
||||
AA2VA1E1E2:
|
||||
SA2VA1E1E2:
|
||||
|
||||
NA2VUE1E2:
|
||||
PA2VUE1E2:
|
||||
AA2VUE1E2:
|
||||
SA2VUE1E2:
|
||||
|
||||
NA2VA1UE1E2:
|
||||
PA2VA1UE1E2:
|
||||
AA2VA1UE1E2:
|
||||
SA2VA1UE1E2:
|
||||
|
||||
NA1VA1:
|
||||
PA1VA1:
|
||||
AA1VA1:
|
||||
SA1VA1:
|
||||
|
||||
NA1VU:
|
||||
PA1VU:
|
||||
AA1VU:
|
||||
SA1VU:
|
||||
|
||||
NA1VE1:
|
||||
PA1VE1:
|
||||
AA1VE1:
|
||||
SA1VE1:
|
||||
|
||||
NA1VE2:
|
||||
PA1VE2:
|
||||
AA1VE2:
|
||||
SA1VE2:
|
||||
|
||||
NA1VA1U:
|
||||
PA1VA1U:
|
||||
AA1VA1U:
|
||||
SA1VA1U:
|
||||
|
||||
NA1VA1E1:
|
||||
PA1VA1E1:
|
||||
AA1VA1E1:
|
||||
SA1VA1E1:
|
||||
|
||||
NA1VA1E2:
|
||||
PA1VA1E2:
|
||||
AA1VA1E2:
|
||||
SA1VA1E2:
|
||||
|
||||
NA1VUE1:
|
||||
PA1VUE1:
|
||||
AA1VUE1:
|
||||
SA1VUE1:
|
||||
|
||||
NA1VUE2:
|
||||
PA1VUE2:
|
||||
AA1VUE2:
|
||||
SA1VUE2:
|
||||
|
||||
NA1VE1E2:
|
||||
PA1VE1E2:
|
||||
AA1VE1E2:
|
||||
SA1VE1E2:
|
||||
|
||||
NA1VA1UE1:
|
||||
PA1VA1UE1:
|
||||
AA1VA1UE1:
|
||||
SA1VA1UE1:
|
||||
|
||||
NA1VA1UE2:
|
||||
PA1VA1UE2:
|
||||
AA1VA1UE2:
|
||||
SA1VA1UE2:
|
||||
|
||||
NA1VA1E1E2:
|
||||
PA1VA1E1E2:
|
||||
AA1VA1E1E2:
|
||||
SA1VA1E1E2:
|
||||
|
||||
NA1VUE1E2:
|
||||
PA1VUE1E2:
|
||||
AA1VUE1E2:
|
||||
SA1VUE1E2:
|
||||
|
||||
NA1VA1UE1E2:
|
||||
PA1VA1UE1E2:
|
||||
AA1VA1UE1E2:
|
||||
SA1VA1UE1E2:
|
||||
|
||||
NEWZE1:
|
||||
PEWZE1:
|
||||
AEWZE1:
|
||||
SEWZE1:
|
||||
|
||||
NEWZE2:
|
||||
PEWZE2:
|
||||
AEWZE2:
|
||||
SEWZE2:
|
||||
|
||||
NEWZE1A1:
|
||||
PEWZE1A1:
|
||||
AEWZE1A1:
|
||||
SEWZE1A1:
|
||||
|
||||
NEWZE2A1:
|
||||
PEWZE2A1:
|
||||
AEWZE2A1:
|
||||
SEWZE2A1:
|
||||
|
||||
NEWZE1U:
|
||||
PEWZE1U:
|
||||
AEWZE1U:
|
||||
SEWZE1U:
|
||||
|
||||
NEWZE2U:
|
||||
PEWZE2U:
|
||||
AEWZE2U:
|
||||
SEWZE2U:
|
||||
|
||||
NEWZE1A1U:
|
||||
PEWZE1A1U:
|
||||
AEWZE1A1U:
|
||||
SEWZE1A1U:
|
||||
|
||||
NEWZE2A1U:
|
||||
PEWZE2A1U:
|
||||
AEWZE2A1U:
|
||||
SEWZE2A1U:
|
||||
|
||||
NEWZE1E2:
|
||||
PEWZE1E2:
|
||||
AEWZE1E2:
|
||||
SEWZE1E2:
|
||||
|
||||
NEWZE1E2A1:
|
||||
PEWZE1E2A1:
|
||||
AEWZE1E2A1:
|
||||
SEWZE1E2A1:
|
||||
|
||||
NEWZE1E2U:
|
||||
PEWZE1E2U:
|
||||
AEWZE1E2U:
|
||||
SEWZE1E2U:
|
||||
|
||||
NEWZE1E2A1U:
|
||||
PEWZE1E2A1U:
|
||||
AEWZE1E2A1U:
|
||||
SEWZE1E2A1U:
|
||||
|
||||
NA1EWZE1:
|
||||
PA1EWZE1:
|
||||
AA1EWZE1:
|
||||
SA1EWZE1:
|
||||
|
||||
NA1EWZE2:
|
||||
PA1EWZE2:
|
||||
AA1EWZE2:
|
||||
SA1EWZE2:
|
||||
|
||||
NA1EWZE1A1:
|
||||
PA1EWZE1A1:
|
||||
AA1EWZE1A1:
|
||||
SA1EWZE1A1:
|
||||
|
||||
NA1EWZE2A1:
|
||||
PA1EWZE2A1:
|
||||
AA1EWZE2A1:
|
||||
SA1EWZE2A1:
|
||||
|
||||
NA1EWZE1U:
|
||||
PA1EWZE1U:
|
||||
AA1EWZE1U:
|
||||
SA1EWZE1U:
|
||||
|
||||
NA1EWZE2U:
|
||||
PA1EWZE2U:
|
||||
AA1EWZE2U:
|
||||
SA1EWZE2U:
|
||||
|
||||
NA1EWZE1A1U:
|
||||
PA1EWZE1A1U:
|
||||
AA1EWZE1A1U:
|
||||
SA1EWZE1A1U:
|
||||
|
||||
NA1EWZE2A1U:
|
||||
PA1EWZE2A1U:
|
||||
AA1EWZE2A1U:
|
||||
SA1EWZE2A1U:
|
||||
|
||||
NA1EWZE1E2:
|
||||
PA1EWZE1E2:
|
||||
AA1EWZE1E2:
|
||||
SA1EWZE1E2:
|
||||
|
||||
NA1EWZE1E2A1:
|
||||
PA1EWZE1E2A1:
|
||||
AA1EWZE1E2A1:
|
||||
SA1EWZE1E2A1:
|
||||
|
||||
NA1EWZE1E2U:
|
||||
PA1EWZE1E2U:
|
||||
AA1EWZE1E2U:
|
||||
SA1EWZE1E2U:
|
||||
|
||||
NA1EWZE1E2A1U:
|
||||
PA1EWZE1E2A1U:
|
||||
AA1EWZE1E2A1U:
|
||||
SA1EWZE1E2A1U:
|
||||
|
||||
NA2EQWZE1U:
|
||||
PA2EQWZE1U:
|
||||
AA2EQWZE1U:
|
||||
SA2EQWZE1U:
|
||||
|
||||
NA2EQWZE2U:
|
||||
PA2EQWZE2U:
|
||||
AA2EQWZE2U:
|
||||
SA2EQWZE2U:
|
||||
|
||||
NA2EQWZE1UA1:
|
||||
PA2EQWZE1UA1:
|
||||
AA2EQWZE1UA1:
|
||||
SA2EQWZE1UA1:
|
||||
|
||||
NA2EQWZE2UA1:
|
||||
PA2EQWZE2UA1:
|
||||
AA2EQWZE2UA1:
|
||||
SA2EQWZE2UA1:
|
||||
|
||||
NA2EQWZE1E2U:
|
||||
PA2EQWZE1E2U:
|
||||
AA2EQWZE1E2U:
|
||||
SA2EQWZE1E2U:
|
||||
|
||||
NA2EQWZE1E2UA1:
|
||||
PA2EQWZE1E2UA1:
|
||||
AA2EQWZE1E2UA1:
|
||||
SA2EQWZE1E2UA1:
|
||||
|
||||
NA2A1QU:
|
||||
PA2A1QU:
|
||||
AA2A1QU:
|
||||
SA2A1QU:
|
||||
|
||||
NA2A1QUA1:
|
||||
PA2A1QUA1:
|
||||
AA2A1QUA1:
|
||||
SA2A1QUA1:
|
||||
|
||||
NA2A1QUE1:
|
||||
PA2A1QUE1:
|
||||
AA2A1QUE1:
|
||||
SA2A1QUE1:
|
||||
|
||||
NA2A1QUE2:
|
||||
PA2A1QUE2:
|
||||
AA2A1QUE2:
|
||||
SA2A1QUE2:
|
||||
|
||||
NA2A1QUA1E1:
|
||||
PA2A1QUA1E1:
|
||||
AA2A1QUA1E1:
|
||||
SA2A1QUA1E1:
|
||||
|
||||
NA2A1QUA1E2:
|
||||
PA2A1QUA1E2:
|
||||
AA2A1QUA1E2:
|
||||
SA2A1QUA1E2:
|
||||
|
||||
NA2A1QUE1E2:
|
||||
PA2A1QUE1E2:
|
||||
AA2A1QUE1E2:
|
||||
SA2A1QUE1E2:
|
||||
|
||||
NA2A1QUA1E1E2:
|
||||
PA2A1QUA1E1E2:
|
||||
AA2A1QUA1E1E2:
|
||||
SA2A1QUA1E1E2:
|
||||
|
||||
NCORECORE:
|
||||
PCORECORE:
|
||||
ACORECORE:
|
||||
SCORECORE:
|
||||
|
||||
NVA2N1:
|
||||
PVA2N1:
|
||||
AVA2N1:
|
||||
SVA2N1:
|
||||
|
||||
NVE1N1:
|
||||
PVE1N1:
|
||||
AVE1N1:
|
||||
SVE1N1:
|
||||
|
||||
NVA1N1:
|
||||
PVA1N1:
|
||||
AVA1N1:
|
||||
SVA1N1:
|
||||
|
||||
NVA2N2:
|
||||
PVA2N2:
|
||||
AVA2N2:
|
||||
SVA2N2:
|
||||
|
||||
NVE1N2:
|
||||
PVE1N2:
|
||||
AVE1N2:
|
||||
SVE1N2:
|
||||
|
||||
NVA1N2:
|
||||
PVA1N2:
|
||||
AVA1N2:
|
||||
SVA1N2:
|
||||
|
||||
NVA2N3:
|
||||
PVA2N3:
|
||||
AVA2N3:
|
||||
SVA2N3:
|
||||
|
||||
NVE1N3:
|
||||
PVE1N3:
|
||||
AVE1N3:
|
||||
SVE1N3:
|
||||
|
||||
NVA1N3:
|
||||
PVA1N3:
|
||||
AVA1N3:
|
||||
SVA1N3:
|
||||
|
||||
NJE1N0:
|
||||
PJE1N0:
|
||||
AJE1N0:
|
||||
SJE1N0:
|
||||
|
||||
NJE1N1:
|
||||
PJE1N1:
|
||||
AJE1N1:
|
||||
SJE1N1:
|
||||
|
||||
NJE1N2:
|
||||
PJE1N2:
|
||||
AJE1N2:
|
||||
SJE1N2:
|
||||
|
||||
NJE1N3:
|
||||
PJE1N3:
|
||||
AJE1N3:
|
||||
SJE1N3:
|
||||
|
||||
NPA2E1N0:
|
||||
PPA2E1N0:
|
||||
APA2E1N0:
|
||||
SPA2E1N0:
|
||||
|
||||
NPA2A1N1:
|
||||
PPA2A1N1:
|
||||
APA2A1N1:
|
||||
SPA2A1N1:
|
||||
|
||||
NPA2E1N2:
|
||||
PPA2E1N2:
|
||||
APA2E1N2:
|
||||
SPA2E1N2:
|
||||
|
||||
NPA2A1N2:
|
||||
PPA2A1N2:
|
||||
APA2A1N2:
|
||||
SPA2A1N2:
|
||||
|
||||
NPA2E1N3:
|
||||
PPA2E1N3:
|
||||
APA2E1N3:
|
||||
SPA2E1N3:
|
||||
|
||||
NPA2A1N3:
|
||||
PPA2A1N3:
|
||||
APA2A1N3:
|
||||
SPA2A1N3:
|
||||
|
||||
NPE1A1N0:
|
||||
PPE1A1N0:
|
||||
APE1A1N0:
|
||||
SPE1A1N0:
|
||||
|
||||
NPE1A1N2:
|
||||
PPE1A1N2:
|
||||
APE1A1N2:
|
||||
SPE1A1N2:
|
||||
|
||||
NPE1A1N3:
|
||||
PPE1A1N3:
|
||||
APE1A1N3:
|
||||
SPE1A1N3:
|
||||
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,757 @@
|
|||
module write_mod
|
||||
implicit none
|
||||
! unit conversion
|
||||
double precision ,parameter :: h2icm = 219474.69d0
|
||||
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 ='20pt_pyram_param.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)
|
||||
! 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*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, '(',f8.1,' icm)')
|
||||
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
|
||||
Loading…
Reference in New Issue