From de65a49694cd0994aa2f8fd6a5dc2b48adf65365 Mon Sep 17 00:00:00 2001 From: jean paul nshuti Date: Wed, 8 Oct 2025 15:23:56 +0200 Subject: [PATCH] first commit of the branch of lmatrix of no3 --- src/model/.giosaveVT1ebr | 0 src/model/JTmod.incl | 38 ++ src/model/adia.f90 | 120 ++++++ src/model/ctrans.f90 | 386 ++++++++++++++++++ src/model/data_transform.f90 | 133 ++++++ src/model/keys.f90 | 84 ++++ src/model/matrix_form.f90 | 374 +++++++++++++++++ src/model/model.f90 | 507 +++++++++++++++++++++++ src/model/nnparams.incl | 43 ++ src/model/surface_mod.f90 | 85 ++++ src/model/weight.f | 50 +++ src/model/write.f | 763 +++++++++++++++++++++++++++++++++++ 12 files changed, 2583 insertions(+) create mode 100644 src/model/.giosaveVT1ebr create mode 100644 src/model/JTmod.incl create mode 100644 src/model/adia.f90 create mode 100644 src/model/ctrans.f90 create mode 100644 src/model/data_transform.f90 create mode 100644 src/model/keys.f90 create mode 100644 src/model/matrix_form.f90 create mode 100644 src/model/model.f90 create mode 100644 src/model/nnparams.incl create mode 100644 src/model/surface_mod.f90 create mode 100644 src/model/weight.f create mode 100644 src/model/write.f diff --git a/src/model/.giosaveVT1ebr b/src/model/.giosaveVT1ebr new file mode 100644 index 0000000..e69de29 diff --git a/src/model/JTmod.incl b/src/model/JTmod.incl new file mode 100644 index 0000000..4092be0 --- /dev/null +++ b/src/model/JTmod.incl @@ -0,0 +1,38 @@ +!*** Relevant parameters for the analytic model +!*** offsets: +!*** offsets(1): morse equilibrium (N-H, Angström) +!*** offsets(2): reference angle (H-N-H) +!*** offsets(3): -- +!*** pat_index: vector giving the position of the +!*** various coordinates (see below) +!*** ppars: polynomial parameters for tmcs +!*** vcfs: coefficients for V expressions. +!*** wzcfs: coefficients for W & Z expressions. +!*** ifc: inverse factorials. + + integer matdim + parameter (matdim=5) ! matrix is (matdim)x(matdim) + + real*8 offsets(2) + integer pat_index(maxnin) + +! NH3 params + parameter (offsets=[2.344419d0,120.d0]) + +!########################################################################## +! coordinate order; the first #I number of coords are given to the +! ANN, where #I is the number of input neurons. The position i in +! pat_index corresponds to a coordinate, the value of pat_index(i) +! signifies its position. +! +! The vector is ordered as follows: +! a,xs,ys,xb,yb,b,rs**2,rb**2,b**2, +! es*eb, es**3, eb**3,es**2*eb, es*eb**2 +! ri**2 := xi**2+yi**2 = ei**2; ei := (xi,yi), i = s,b +! +! parts not supposed to be read by ANN are marked by ';' for your +! convenience. +!########################################################################## +! a,rs**2,rb**2,es*eb,es**3,eb**3,es**2*eb,es*eb**2,b**2 #I=9 (6D) + parameter (pat_index=[1,2,3,4,5,6,7,8,9,10,11,12,13,14]) +!************************************************************************** diff --git a/src/model/adia.f90 b/src/model/adia.f90 new file mode 100644 index 0000000..78a0c91 --- /dev/null +++ b/src/model/adia.f90 @@ -0,0 +1,120 @@ + module adia_mod + implicit none + contains + +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +! % SUBROUTINE ADIA(N,P,NPAR,ymod,v,u,SKIP) +! % +! % determines the adiabatic energies by diagonalizing diabatic matrix. +! % The Eingenvalues are sorted according to the best fitting ordering +! % of the CI vectors. +! % +! % ATTENTION: The interface has changed. To sort by the ci's, +! % the datavalue of the current points are given +! % +! % input variables: +! % n: number of point (int) +! % p: parameter evector(double[npar]) +! % npar: number of parameters (int) +! % skip: .false. if everything should be done +! % +! % output variables: +! % ymod: firtst nstat energies and than nci*ndiab ci's (double[ntot]) +! % v: eigenvalues (double[ndiab]) +! % u: eigenvectors (double[ndiab,ndiab]) +!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + subroutine adia(n,p,npar,ymod,vx,u,skip) + use dim_parameter,only: ndiab,nstat,ntot,nci,pst + use data_module,only: q_m,x1_m,x2_m,y_m + use diab_mod, only:diab + use data_matrix + !use dipole, only: diab + implicit none + + integer i,j !running indices + + integer iref !getting correction or refference + + double precision e(ndiab,ndiab) !full diabatic matrix + double precision mx(ndiab,ndiab) + double precision my(ndiab,ndiab) + double precision vxs,vys,vxb,vyb + + integer n !current point + integer npar !number of parameters + double precision p(npar) !parameters + + double precision u(ndiab,ndiab),ut(ndiab,ndiab) !ci-vectors + double precision ymod(ntot) !fitted data + double precision vx(ndiab),vy(nstat) !eigen values + double precision,allocatable,dimension(:,:):: mat + logical skip,dbg + parameter (dbg=.false.) + double precision,dimension(2,2):: T,TT,TX,TY +! lapack variables + integer,parameter :: lwork = 1000 + double precision work(lwork) + integer info + integer TYPES, BLK ! TYPE OF THE CALCULATION + ! variabke for dgemm + + double precision,dimension(ndiab,ndiab):: ex,ey,ez + double precision:: alpha + integer:: lda,ldb,beta,ldc + double precision,dimension(ndiab,ndiab):: temp1,temp2 + call diab(ex,ey,ez,n,x1_m(:,n),x2_m(:,n),p) + + +! init eigenvector matrix + TYPES = int(p(pst(1,33))) + + BLK = int(p(pst(1,33)+1)) ! BLOCK IF TYPE IS 3 + u = 0.d0 + vx=0.0d0 + skip=.false. + ymod=0.0d0 + if (TYPES .eq.1 ) then + ! Trace of the potential + call trace_mat(ex,ey,ymod) + else if (TYPES .eq.2) then + ! Eigenvalue decomposition of the potential + call Eigen(ex,ey,ymod) + + else if (TYPES .eq.3) then + CALL BLOCK_DIAB(ex,ey,ymod,BLK) + else if (TYPES .EQ.4) then + call Full_diab_upper(ex,ey,ymod) + else if (TYPES .eq.5) then + call Transformation_mat(ex,vx,ymod) + ymod=0.0d0 + else if (TYPES .eq.6) then + ! transform the lz + call one_dia_upper(ez,ymod) + else + write(*,*) "Error in TYPE of calculation here",TYPES + + stop + end if + if (dbg) then + do i=1,ndiab + write(*,'(5f14.6)') (ex(i,j),j=1,ndiab) + enddo + write(*,*)"" + endif + + end subroutine + subroutine matrix_mult(C,A,B,N) + implicit none + integer:: n,i,j,k + double precision,dimension(n,n):: A,B,C + do i = 1, n ! Rows of C + do j = 1, n ! Columns of C + C(i,j) = 0.0 ! Initialize element + do k = 1, n ! Dot product + C(i,j) = C(i,j) + A(i,k) * B(k,j) + end do + end do + end do + end subroutine + + end module adia_mod diff --git a/src/model/ctrans.f90 b/src/model/ctrans.f90 new file mode 100644 index 0000000..d50bdb9 --- /dev/null +++ b/src/model/ctrans.f90 @@ -0,0 +1,386 @@ +module ctrans_mod + use dim_parameter, only: qn + contains +!! subroutine ctrans + + subroutine ctrans(q,x1,x2) + implicit none + include 'nnparams.incl' + include 'JTmod.incl' + double precision,intent(in):: q(qn) + double precision,intent(out):: x1(qn),x2(qn) + double precision:: cart(3,4),qint(maxnin) + integer i + !cart(:,1)=0.0d0 + + !cart(1:3,2:4) = reshape([ q(4:12) ], shape(cart(1:3,2:4))) + cart(1,1)=q(1) + cart(2,1)=q(2) + cart(3,1)=q(3) + cart(1,2)=q(4) + cart(2,2)=q(5) + cart(3,2)=q(6) + cart(1,3)=q(7) + cart(2,3)=q(8) + cart(3,3)=q(9) + cart(1,4)=q(10) + cart(2,4)=q(11) + cart(3,4)=q(12) + call cart2int(cart,qint) + do i=1,qn + if (abs(qint(i)) .lt. 1.0d-5) qint(i) =0.0d0 + enddo + x1(1:qn)=qint(1:qn) + !x1(2)=0.0d0 + x1(5)=-x1(5) + x1(3)=-x1(3) + !x1(6)=0.0d0 + x2(1:qn)=0.0d0 !qint(1:qn) + end subroutine ctrans + + + subroutine cart2int(cart,qint) + implicit none +! This version merges both coordinate transformation routines into +! one. JTmod's sscales(2:3) are ignored. +! This is the first version to be compatible with one of my proper 6D fits +! Time-stamp: <2024-10-22 13:52:59 dwilliams> + +! Input (cartesian, in Angström) +! cart(:,1): N +! cart(:,1+i): Hi +! Output +! qint(i): order defined in JTmod. +! Internal Variables +! no(1:3): NO distances 1-3 +! pat_in: temporary coordinates +! axis: main axis of NO3 + include 'nnparams.incl' + include 'JTmod.incl' + + + real*8 cart(3,4),qint(maxnin) + + real*8 no(3), r1, r2, r3 + real*8 v1(3), v2(3), v3(3) + real*8 n1(3), n2(3), n3(3), tr(3) + real*8 ortho(3) + real*8 pat_in(maxnin) + logical ignore_umbrella,dbg_umbrella + logical dbg_distances + +!.. Debugging parameters +!.. set umbrella to 0 + parameter (ignore_umbrella=.false.) +! parameter (ignore_umbrella=.true.) + +!.. break if umbrella is not 0 + parameter (dbg_umbrella=.false.) +! parameter (dbg_umbrella=.true.) + +!.. break for tiny distances + parameter (dbg_distances=.false.) +! parameter (dbg_distances=.true.) + + integer k + +!.. get N-O vectors and distances: + do k=1,3 + v1(k)=cart(k,2)-cart(k,1) + v2(k)=cart(k,3)-cart(k,1) + v3(k)=cart(k,4)-cart(k,1) + enddo + no(1)=norm(v1,3) + no(2)=norm(v2,3) + no(3)=norm(v3,3) + +!.. temporarily store displacements + do k=1,3 + pat_in(k)=no(k)-offsets(1) + enddo + + do k=1,3 + v1(k)=v1(k)/no(1) + v2(k)=v2(k)/no(2) + v3(k)=v3(k)/no(3) + enddo + +!.. compute three normal vectors for the ONO planes: + call xprod(n1,v1,v2) + call xprod(n2,v2,v3) + call xprod(n3,v3,v1) + + do k=1,3 + tr(k)=(n1(k)+n2(k)+n3(k))/3.d0 + enddo + r1=norm(tr,3) + do k=1,3 + tr(k)=tr(k)/r1 + enddo + +! rotate trisector + call rot_trisec(tr,v1,v2,v3) + +!.. determine trisector angle: + if (ignore_umbrella) then + pat_in(7)=0.0d0 + else + pat_in(7)=pi/2.0d0 - acos(scalar(v1,tr,3)) + pat_in(7)=sign(pat_in(7),cart(1,2)) + endif + +!.. molecule now lies in yz plane, compute projected ONO angles: + v1(1)=0.d0 + v2(1)=0.d0 + v3(1)=0.d0 + r1=norm(v1,3) + r2=norm(v2,3) + r3=norm(v3,3) + do k=2,3 + v1(k)=v1(k)/r1 + v2(k)=v2(k)/r2 + v3(k)=v3(k)/r3 + enddo + +! make orthogonal vector to v3 + ortho(1)=0.0d0 + ortho(2)=v3(3) + ortho(3)=-v3(2) + +!.. projected ONO angles in radians + pat_in(4)=get_ang(v2,v3,ortho) + pat_in(5)=get_ang(v1,v3,ortho) + + pat_in(6)=dabs(pat_in(5)-pat_in(4)) + +!.. account for rotational order of atoms + if (pat_in(4).le.pat_in(5)) then + pat_in(5)=2*pi-pat_in(4)-pat_in(6) + else + pat_in(4)=2*pi-pat_in(5)-pat_in(6) + endif + + pat_in(4)=rad2deg*pat_in(4)-offsets(2) + pat_in(5)=rad2deg*pat_in(5)-offsets(2) + pat_in(6)=rad2deg*pat_in(6)-offsets(2) + pat_in(7)=rad2deg*pat_in(7) + + call genANN_ctrans(pat_in) + + qint(:)=pat_in(:) + + contains + +!------------------------------------------------------------------- +! compute vector product n1 of vectors v1 x v2 + subroutine xprod(n1,v1,v2) + implicit none + + real*8 n1(3), v1(3), v2(3) + + n1(1) = v1(2)*v2(3) - v1(3)*v2(2) + n1(2) = v1(3)*v2(1) - v1(1)*v2(3) + n1(3) = v1(1)*v2(2) - v1(2)*v2(1) + + end subroutine + + +!------------------------------------------------------------------- +! compute scalar product of vectors v1 and v2: + real*8 function scalar(v1,v2,n) + implicit none + integer i, n + real*8 v1(*), v2(*) + + scalar=0.d0 + do i=1,n + scalar=scalar+v1(i)*v2(i) + enddo + + end function + +!------------------------------------------------------------------- +! compute norm of vector: + real*8 function norm(x,n) + implicit none + integer i, n + real*8 x(*) + + norm=0.d0 + do i=1,n + norm=norm+x(i)**2 + enddo + norm=sqrt(norm) + + end function + +!------------------------------------------------------------------- + + subroutine rot_trisec(tr,v1,v2,v3) + implicit none + + real*8 tr(3),v1(3),v2(3),v3(3) + + real*8 vrot(3) + real*8 rot_ax(3) + real*8 cos_phi,sin_phi + +! evaluate cos(-phi) and sin(-phi), where phi is the angle between +! tr and (1,0,0) + cos_phi=tr(1) + sin_phi=dsqrt(tr(2)**2+tr(3)**2) + + if (sin_phi.lt.1.0d-12) then + return + endif + +! determine rotational axis + rot_ax(1) = 0.0d0 + rot_ax(2) = tr(3) + rot_ax(3) = -tr(2) +! normalize + rot_ax=rot_ax/sin_phi + +! now the rotation can be done using Rodrigues' rotation formula +! v'=v*cos(p) + (k x v)sin(p) + k (k*v) (1-cos(p)) +! for v=tr k*v vanishes by construction: + +! check that the rotation does what it should + call rodrigues(vrot,tr,rot_ax,cos_phi,sin_phi) + if (dsqrt(vrot(2)**2+vrot(3)**2).gt.1.0d-12) then + write(6,*) "ERROR: BROKEN TRISECTOR" + stop + endif + tr=vrot + + call rodrigues(vrot,v1,rot_ax,cos_phi,sin_phi) + v1=vrot + call rodrigues(vrot,v2,rot_ax,cos_phi,sin_phi) + v2=vrot + call rodrigues(vrot,v3,rot_ax,cos_phi,sin_phi) + v3=vrot + + + end subroutine + +!------------------------------------------------------------------- + + subroutine rodrigues(vrot,v,axis,cos_phi,sin_phi) + implicit none + real*8 vrot(3),v(3),axis(3) + real*8 cos_phi,sin_phi + + real*8 ortho(3) + + call xprod(ortho,axis,v) + vrot = v*cos_phi + ortho*sin_phi+axis*scalar(axis,v,3)*(1-cos_phi) + + end subroutine + +!------------------------------------------------------------------- + + real*8 function get_ang(v,xaxis,yaxis) + implicit none +! get normalized [0:2pi) angle from vectors in the yz plane + real*8 v(3),xaxis(3),yaxis(3) + + real*8 phi + + real*8 pi + parameter (pi=3.141592653589793d0) + + phi=atan2(scalar(yaxis,v,3),scalar(xaxis,v,3)) + if (phi.lt.0.0d0) then + phi=2*pi+phi + endif + get_ang=phi + + end function + + end subroutine cart2int + subroutine genANN_ctrans(pat_in) + implicit none + + include 'nnparams.incl' + include 'JTmod.incl' + + real*8 pat_in(maxnin) + + real*8 raw_in(maxnin),off_in(maxnin),ptrans_in(7) + real*8 r0 + real*8 a,b,xs,ys,xb,yb + + integer k + + off_in(1:7)=pat_in(1:7) + r0=offsets(1) + +! transform primitives +! recover raw distances from offset coords + do k=1,3 + raw_in(k)=off_in(k)+offsets(1) + enddo + + do k=1,3 + ptrans_in(k)=off_in(k) + enddo + +! rescale ONO angles + ptrans_in(4)=deg2rad*off_in(4) + ptrans_in(5)=deg2rad*off_in(5) + ptrans_in(6)=deg2rad*off_in(6) +! rescale umbrella + ptrans_in(7)=off_in(7)*deg2rad + +! compute symmetry coordinates + +! A (breathing) + a=(ptrans_in(1)+ptrans_in(2)+ptrans_in(3))/dsqrt(3.0d0) +! ES + call prim2emode(ptrans_in(1:3),xs,ys) +! EB + call prim2emode(ptrans_in(4:6),xb,yb) +! B (umbrella) + b=ptrans_in(7) + +! overwrite input with output + + pat_in(pat_index(1))=a ! 1 + pat_in(pat_index(2))=xs + pat_in(pat_index(3))=ys + pat_in(pat_index(4))=xb + pat_in(pat_index(5))=yb + pat_in(pat_index(6))=b +! totally symmetric monomials + pat_in(pat_index(7))=xs**2 + ys**2 ! 2 + pat_in(pat_index(8))=xb**2 + yb**2 ! 3 + pat_in(pat_index(9))=b**2 ! 9 + pat_in(pat_index(10))=xs*xb+ys*yb ! 4 +! S^3, B^3 + pat_in(pat_index(11))=xs*(xs**2-3*ys**2) ! 5 + pat_in(pat_index(12))=xb*(xb**2-3*yb**2) ! 6 +! S^2 B, S B^2 + pat_in(pat_index(13))=xb*(xs**2-ys**2) - 2*yb*xs*ys ! 7 + pat_in(pat_index(14))=xs*(xb**2-yb**2) - 2*ys*xb*yb ! 8 + + do k=11,14 + pat_in(pat_index(k))=tanh(0.1d0*pat_in(pat_index(k)))*10.0d0 + enddo + + + + + end subroutine + subroutine prim2emode(prim,ex,ey) + implicit none +! Takes a 2D-vector prim and returns the degenerate modes x and y +! following our standard conventions. + + real*8 prim(3),ex,ey + + ex=(2.0d0*prim(1)-prim(2)-prim(3))/dsqrt(6.0d0) + ey=(prim(2)-prim(3))/dsqrt(2.0d0) + + end +end module ctrans_mod + diff --git a/src/model/data_transform.f90 b/src/model/data_transform.f90 new file mode 100644 index 0000000..411a09a --- /dev/null +++ b/src/model/data_transform.f90 @@ -0,0 +1,133 @@ +! 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 + +!---------------------------------------------------------------------------------------------------- +! (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 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 + +!---------------------------------------------------------------------------------------------------- +! ' Debye)') + 42 format('#',10x,'No. of Points: ',i10) + 43 format('#',10x,'Total weighted RMS: ',g16.8, '(',f8.1,' icm)') + 44 format('#',10x,'Sum of point weights: ',f16.8) + 45 format('#',10x,'Total Energie RMS: ',g16.8, '(',f8.1,' icm)') + + 48 format('#',10x,'Red. Energie RMS: ',g16.8,'(',f8.1,' icm)') + 51 format('#') + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotfiles(x,y,wt,ymod) + use dim_parameter,only: ndata,sets,qn,ntot,numdatpt,plot_coord + implicit none +! IN: variables + double precision x(qn,numdatpt),y(ntot,numdatpt) + double precision wt(ntot,numdatpt), ymod(ntot,numdatpt) +! INTERNAL: variables + integer sstart,ssend,set,id_plot + +! Initialize position pointer + ssend=0 +! loop over datasets and print the plotfiles + do set=1 ,sets + if(ndata(set).eq.0) cycle + id_plot=50+set + call open_plotfile(id_plot,set) + write(id_plot,'(A)') '# -*- truncate-lines: t -*-' +! get start and end point of each set + sstart=ssend+1 + ssend=ssend+ndata(set) + if (plot_coord(set).eq.0) then + call print_plotwalk(x(:,sstart:ssend),y(:,sstart:ssend), + > wt(:,sstart:ssend),ymod(:,sstart:ssend), + > ndata(set),id_plot,set) + else + call print_plotcoord(plot_coord(set), + > x(:,sstart:ssend),y(:,sstart:ssend), + > wt(:,sstart:ssend),ymod(:,sstart:ssend), + > ndata(set),id_plot,set) + endif + close(id_plot) + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotwalk(x,y,wt,ymod,npt,id_plot,set) + use dim_parameter,only: qn,ntot + use io_parameters,only: llen + implicit none +! IN: variables + integer id_plot,npt,set + double precision x(qn,npt),y(ntot,npt),ymod(ntot,npt),wt(ntot,npt) +! INTERNAL: variables + double precision xdiff(qn),walktime + double precision walknorm +! loop control + integer i,j + + character(len=llen) fmt + j=ntot-1 + + call print_plotheader(id_plot,0,npt,set) + + call getwalknorm(x,walknorm,npt) + walktime = 0.d0 + do i=1,npt + if(i.gt.1) then + xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1) + walktime = walktime + dsqrt(sum(xdiff(1:qn)**2))/walknorm + endif + write(id_plot,"(ES16.8,*(3(ES16.8),:))") + > walktime ,ymod(:,i),y(:,i),(wt(:,i)) + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotcoord(coord,x,y,wt,ymod,npt,id_plot,set) + use dim_parameter,only: qn,ntot + use io_parameters,only: llen + implicit none +! IN: variables + integer, intent(in) :: id_plot,npt,set,coord + double precision, intent(in) :: x(qn,npt),y(ntot,npt) + double precision, intent(in) :: ymod(ntot,npt),wt(ntot,npt) +! loop control + integer i + + call print_plotheader(id_plot,coord,npt,set) + do i=1,npt +! write(id_plot,"(ES16.8,*(3(ES16.8),:))") +! > x(coord,i), ymod(:,i),y(:,i),(wt(:,i)) + write(id_plot,"(2ES16.8,*(3(ES16.8),:))") + > x(coord,i), x(coord+1,i),y(:,i) + enddo + + end subroutine + +!---------------------------------------------------------------------------------------------------- + subroutine print_plotheader(id_plot,coord,npt,set) + use dim_parameter,only: qn,ntot + use io_parameters,only: llen + implicit none + integer, intent(in) :: id_plot,npt,set,coord + + character(len=llen) fmt + + write(id_plot,'("#SET: ",i5)') set + write(id_plot,'("#OUTPUT VALUES",i4)') ntot + write(id_plot,'("#DATA POINTS: ",i4)') npt + if (coord.le.0) then + write(id_plot,'("#t(x) = WALK")') + else + write(id_plot,'("#t(x) = x(",I0,")")') coord + endif + write(id_plot,'("#UNIT: hartree")') + write(id_plot,'()') + write(id_plot,'("#",A15)',advance='no') "t(x)" + write(fmt,'("(3(7X,A9,",I3,"(16x)))")') ntot-1 + write(id_plot,fmt) 'ymod(p,x)','y(x) ','wt(x) ' + + + end subroutine + +!---------------------------------------------------------------------------------------------------- +! 'ERROR: No rule for Outputfile naming for MAXIT >= 1000' + stop + endif + + open (id_out,file=outname) + + end subroutine + +!---------------------------------------------------------------------------------------------------- +! =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 + +!---------------------------------------------------------------------------------------------------- +! 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 +!---------------------------------------------------------------------------------------------------- +! wt(k), x(1:qn_read) + enddo + enddo + + end subroutine + + + + end module write_mod