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

Compare commits

base: jnshuti:210f507f58616f0d30f17b2600c46daf1b051c08
Branches Tags
jnshuti:main
jnshuti:model_muxy_nh3
jnshuti:model_lmat_no3
..
compare: jnshuti:464dca1bbb7c488dfa727068a8df0fb0fa48e25c
Branches Tags
jnshuti:model_muxy_nh3
jnshuti:main
jnshuti:model_lmat_no3

2 Commits
210f507f58 ... 464dca1bbb

Author SHA1 Message Date
jean paul nshuti 464dca1bbb remove PES 2025-10-08 15:52:19 +02:00
jean paul nshuti 733a5d4172 first commit of the model of nh3+ 2025-10-08 15:50:50 +02:00
15 changed files with 3268 additions and 0 deletions
Show Stats Expand all files Collapse all files

0
src/model/.giosaveVT1ebr Normal file
View File

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

@ -0,0 +1,38 @@
!*** Relevant parameters for the analytic model
!*** offsets:
!*** offsets(1): morse equilibrium (N-H, Angström)
!*** offsets(2): reference angle (H-N-H)
!*** offsets(3): --
!*** pat_index: vector giving the position of the
!*** various coordinates (see below)
!*** ppars: polynomial parameters for tmcs
!*** vcfs: coefficients for V expressions.
!*** wzcfs: coefficients for W & Z expressions.
!*** ifc: inverse factorials.
integer matdim
parameter (matdim=5) ! matrix is (matdim)x(matdim)
real*8 offsets(2)
integer pat_index(maxnin)
! NH3 params
parameter (offsets=[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])
!**************************************************************************

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

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

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

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

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

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

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

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

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

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

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

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

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

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

104
src/model/old_keys.f90 Normal file
View File

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

357
src/model/pyramidal_model.f90 Normal file
View File

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

110
src/model/temp Normal file
View File

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

540
src/model/temp.keys Normal file
View File

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

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

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

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

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