diff --git a/src/model/.giosaveVT1ebr b/src/model/.giosaveVT1ebr
deleted file mode 100644
index e69de29..0000000
diff --git a/src/model/JTmod.incl b/src/model/JTmod.incl
deleted file mode 100644
index 4092be0..0000000
--- a/src/model/JTmod.incl
+++ /dev/null
@@ -1,38 +0,0 @@
-!***   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
deleted file mode 100644
index 78a0c91..0000000
--- a/src/model/adia.f90
+++ /dev/null
@@ -1,120 +0,0 @@
-      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
deleted file mode 100644
index d50bdb9..0000000
--- a/src/model/ctrans.f90
+++ /dev/null
@@ -1,386 +0,0 @@
-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
deleted file mode 100644
index 411a09a..0000000
--- a/src/model/data_transform.f90
+++ /dev/null
@@ -1,133 +0,0 @@
-!     <subroutine for manipulating the input Data before the Fit
-      subroutine data_transform(q,x1,x2,y,wt,p,npar,p_act)
-      use accuracy_constants, only: dp,idp
-      use dim_parameter,only : nstat,pst,ntot,qn,numdatpt,ndiab,ndata,sets
-      use ctrans_mod, only: ctrans
-     
-      use surface_mod, only: eval_surface
-      use data_matrix
-!      use david_ctrans_mod, only: ctrans_d
-      implicit none
-!     IN: variables
-      integer(idp),intent(in) :: npar 
-      Real(dp),intent(in) :: q(qn,numdatpt)
-      Real(dp),intent(in) :: p(npar)
-      integer(idp),intent(in) :: p_act(npar)
-
-      ! INOUT: variables 
-      Real(dp),intent(inout) :: y(ntot,numdatpt)
-      Real(dp),intent(inout) :: wt(ntot,numdatpt)
-
-      ! OUT: vairables 
-      Real(dp), intent(out) :: x1(qn,numdatpt),x2(qn,numdatpt)
-      ! internal variables 
-
-      
-      Real(dp),dimension(ndiab,ndiab)::mat_x,mat_y,mat_z,U,V
-  
-      Real(dp),dimension(nstat) :: E
-      integer(idp) pt,i,j,k,l, TYPES, BLK  ! types is for the type of calculation 
-      ! blk is for which block to fit 
-      logical,parameter::  dbg = .false.
-    
-
-       
-      if (pst(2,33) .ne. 2) then  
-
-            write(*,*) "Error in Paramater Keys, TYPE_CAL should be 2 parameter", pst(2,33)
-            stop
-      end if 
-
-      TYPES = int(p(pst(1,33)))! TYPE OF THE CALCULATION
-      BLK= int(p(pst(1,33)+1))! BLOCK IF TYPE IS 3
-      write(*,*) "TYPE of calculation:",TYPES
-
-      pt=1
-
-      do i=1,sets ! loop over the number of sets
-        do j=1,ndata(i) ! loop over the nbr of  points in each sets
-            ! remember to increment pt at the end of the loop 
-          call ctrans(q(1:qn,pt),x1(:,pt),x2(:,pt)) ! transform the coordinate
-
-          ! get the reference U matrix
-
-          !if (j .eq. 3) then 
-          !  call eval_surface(E,V,U_ref,q(1:qn,pt))
-          !  call transform_U(U_ref)
-          !endif 
-
-         
-
-      
-      !do pt=1,numdatpt
-         !call ctrans(q(1:qn,pt),x1(:,pt),x2(:,pt))!  ctrans the dipole cooordinate.
-         write(7,'(I3,*(E17.8))') pt,x1(:,pt)
-
-         call eval_surface(E,V,U,q(1:qn,pt)) 
-         
-         
-       ! Transform U mmatrix
-         call transform_U(U) ! Transform the U matrix
-         ! write U matrix on f16
-          if (dbg) then 
-           !write(7,*) "U matrix at point", pt
-           do k=1,ndiab
-             write(50+i,'(2E17.8,5X,5E17.8)')x1(2:3,pt),(U(k,l),l=1,ndiab)
-           enddo
-           write(50+i,*) ""
-          endif
-         !call overlap(U_ref,U)
-       call Y2mat(y(1:ntot,pt),mat_x,mat_y,mat_z)
-
-       if (TYPES .eq.1 ) then 
-          ! Trace of the potential
-          call trace_mat(mat_x,mat_y,y(1:ntot,pt))
-       else if (TYPES .eq.2) then 
-         ! Eigenvalue decomposition of the potential
-         call Eigen(mat_x,mat_y,y(1:ntot,pt))
-       else if (TYPES .eq.3) then 
-         ! Adiabatic transformation of the potential
-         call adiabatic_transform(mat_x,U)
-         call adiabatic_transform(mat_y,U)
-         call block_diab(mat_x,mat_y,y(1:ntot,pt),BLK)
-         
-       else if (TYPES .eq.4) then
-          ! Write the full upper diabatic matrix 
-          call adiabatic_transform(mat_x,U)
-          call adiabatic_transform(mat_y,U)
-          ! and write the full diabatic matrix to y
-          ! This is the full diabatic matrix
-          call Full_diab_upper(mat_x,mat_y,y(1:ntot,pt))
-        else if (TYPES .eq.5) then
-          !call adiabatic_transform(mat_x,U)
-          !call adiabatic_transform(mat_y,U)
-          call Transformation_mat(U,E,y(1:ntot,pt))
-        else if (TYPES .eq.6) then
-          ! Just do the adiabatic transformation and write the matrix
-          ! transform the lz 
-          call adiabatic_transform(mat_z,U)
-          call one_dia_upper(mat_z,y(1:ntot,pt))
-       else 
-         write(*,*) "Error in TYPE of calculationss",TYPES
-         write(*,*) "the value:,", p(pst(1,33))
-         stop
-       end if
-       pt=pt+1
-      enddo ! j
-      write(34,*) "#---- End of set ", i
-
-      write(7,*) "#---- End of set ", i
-      
-    enddo ! i
-
-      !enddo
-
-      call weight(wt,y)
-      end subroutine
-      
-
-
-      
-     
-
-      
diff --git a/src/model/keys.f90 b/src/model/keys.f90
deleted file mode 100644
index 4cffdd9..0000000
--- a/src/model/keys.f90
+++ /dev/null
@@ -1,84 +0,0 @@
-module keys_mod
-    implicit none 
-contains
-!program gen_key
-!    implicit none
-!    call init_keys()
-!end program gen_key
-    subroutine init_keys
-        use io_parameters, only: key 
-        character(len=1) prefix(4)
-        parameter (prefix=['N','P','A','S'])
-        !character (len=20) key(4,25)
-        integer,parameter:: np=33
-        character(len=16) parname(np)
-        integer i,j
-! Defining keys for potential 
-      ! the electronic state of NO3 A2' E" E'
-    ! Naming convention 
-    ! the keys for Lx and Ly
-    ! the coupling between A2' and A2"
-    parname(1)='LXYVA2O1'
-    parname(2)='LXYVE1O1'
-    parname(3)='LXYVE2O1'
-    parname(4)='LXYVA2O2'
-    parname(5)='LXYVE1O2'
-    parname(6)='LXYVE2O2'
-
-    ! W & Z of E1
-    parname(7)='LXYWZE1O0'
-    parname(8)='LXYWZE1O1'
-    parname(9)='LXYWZE1O2'
-    parname(10)='LXYWZE2O0'
-    parname(11)='LXYWZE2O1'
-    parname(12)='LXYWZE2O2'
-   
-
-    ! WW and Z Pseudo between E1 and E2
-    ! p STANDS FOR PSEUDO JAHN-TELLER
-    parname(13)='LXYPE1E2O0'
-    parname(14)='LXYPE1E2O1'
-    parname(15)='LXYPE1E2O2'
-    ! no order 3
-
-    ! PSEUDO A2 & E1
-    parname(16)='LXYPA2E1O0'
-    parname(17)='LXYPA2E1O1'
-    parname(18)='LXYPA2E1O2'
-
-    ! Pseudo JAHN-TELLER BETWEEN A2 AND E1
-
-    parname(19)='LXYPA2E2O0'
-    parname(20)='LXYPA2E2O1'
-    parname(21)='LXYPA2E2O2'
-    
-
-    ! keys for lz 
-
-    parname(22)='LZWZE1O1'
-    parname(23)='LZWZE1O2'
-    parname(24)='LZWZE2O1'
-    parname(25)='LZWZE2O2'
-    parname(26)='LZPE1E2O0'
-    parname(27)='LZPE1E2O1'
-    parname(28)='LZPE1E2O2'
-    parname(29)='LZPA2E1O1'
-    parname(30)='LZPA2E1O2'
-    parname(31)='LZPA2E2O1'
-    parname(32)='LZPA2E2O2'
-
-
-    parname(33)='TYPE_CAL'! TYPE OF THE CALCULATION WHETHER IT IS THE TRACE OR SOMETHING ELSE
-
-        do i=1,np
-            do j=1,4
-                key(j, i)=prefix(j)//trim(parname(i))//':'  
-             write(8,*) key(j,i)
-            enddo 
-            write(8,*) ''
-           
-        enddo
-       
-        end subroutine
-        
-end module keys_mod
diff --git a/src/model/matrix_form.f90 b/src/model/matrix_form.f90
deleted file mode 100644
index 18fddf1..0000000
--- a/src/model/matrix_form.f90
+++ /dev/null
@@ -1,374 +0,0 @@
- module data_matrix
- use dim_parameter, only:ndiab,nstat,ntot,pst
-! use surface_mod, only: eval_surface
- contains
-      ! subroutine trace 
-
-      subroutine trace_mat(mx,my,y)
-            IMPLICIT NONE 
-            integer::i 
-            double precision,intent(inout):: y(:)
-            double precision, intent(in):: mx(:,:),my(:,:)
-            y=0.0d0
-            !y(1)=mx(4,4)+mx(5,5)
-
-            do i=1,ndiab 
-                  y(1)=y(1)+mx(i,i)
-                  y(2)=y(2)+my(i,i)
-            enddo 
-
-      END SUBROUTINE trace_mat
-            !! subroutine Ydata to matrix
-
-      subroutine Y2mat(Y,Mx,My,mz)
-            IMPLICIT NONE
-            integer:: ii,i,j
-            double precision, intent(in):: y(:)
-            double precision,dimension(ndiab,ndiab),intent(out):: mx, my,mz 
-
-            !if (ndiab .ne. 4 ) then
-            !write(*,*) " NDIAB should be equal to 4",NDIAB
-            !write(*,*) "CHECK DATA_TRANSFORM TO MAKE IT ADAPTABLE"
-            !stop 
-            !endif 
-            ii=1
-            do i=1,ndiab 
-            do j=1,i
-            !      !mx
-
-                  mx(i,j)=y(ii)
-            !      ! My
-                  my(i,j)=y( (ntot/3)+ii)
-                   ! remember to adjust here I added the energy
-                  mz(i,j)= y(2*(ntot/3)+ ii )
-            !      
-                  ii=ii+1
-            enddo  
-            enddo 
-            call copy_2_upper(mx)
-            call copy_2_upper(my)
-            call copy_2_upper(mz)
-      end subroutine
-
-      subroutine Full_diab_upper(mx,my,y)
-            implicit none
-            double precision,intent(inout) :: y(:)
-            double precision, intent(in) :: mx(ndiab,ndiab), my(ndiab,ndiab)
-            integer i,j,ii
-            ii=1
-            y=0.0d0
-            
-            do i=1,ndiab
-                  do j=i,ndiab
-                        ! mx  
-                        y(ii) = mx(i,j)
-                        ! my
-                        y((ntot/2)+ii) = my(i,j)
-                        ! increment the index
-                        ii=ii+1
-                  enddo
-            enddo
-      end subroutine Full_diab_upper
-
-      subroutine one_dia_upper(m,y)
-            implicit none
-            double precision,intent(inout) :: y(:)
-            double precision, intent(in) :: m(ndiab,ndiab)
-            integer i,j,ii
-            ii=1
-            y=0.0d0
-            
-            do i=1,ndiab
-                  do j=i,ndiab
-                        ! mx  
-                        y(ii) = m(i,j)
-                        ! increment the index
-                        ii=ii+1
-                  enddo
-            enddo
-      end subroutine one_dia_upper
-
-
-      Subroutine adiabatic_transform(mx,U)
-            implicit none
-            double precision, intent(inout) :: mx(ndiab,ndiab)
-            double precision, dimension(:,:), intent(inout) :: U
-            double precision, dimension(ndiab,ndiab) :: temp1, temp2
-            integer i, j    
-            
-            !call transform_U(U) ! Transform the U matrix
-
-            ! Transform mx and my to adiabatic basis
-            temp1 = matmul(mx, transpose(U))
-            mx = matmul(U, temp1)         
-            !temp2 = matmul(my, transpose(U))
-            !my = matmul(U, temp2)
-
-      end subroutine adiabatic_transform
-
-       ! the eigenvalue of the dipole 
-
-      SUBROUTINE Eigen(mx,my,Yres)
-            implicit none
-            double precision,dimension(:,:),intent(inout) :: mx,my
-            double precision,dimension(:),intent(out) :: Yres
-            double precision,dimension(ndiab) :: vx,vy
-            double precision,dimension(size(mx,1),size(my,2)) :: temp
-            ! create a temorary matrix fo the eigenvctors
-
-            double precision, allocatable :: mux(:,:), muy(:,:)
-            
-            ! Lapak parameters
-            integer :: n,info,i 
-            integer,parameter :: lwork = 100
-            double precision :: work(lwork)  
-            ! temporary 
-            double precision:: max_row
-            Yres = 0.0d0
-            Allocate(mux,source=mx)
-            call DSYEV('V', 'U', size(mx,1), mux, size(mx,1), vx, work, lwork, info)
-            mx=mux
-            if (info /= 0) then
-            write(*,*) "Error in Eigenvalue decomposition of mx info = ", info
-            stop
-            end if
-            deallocate(mux)
-            Allocate(muy,source=my)
-            call DSYEV('V', 'U', size(my,1), muy, size(my,1), vy, work, lwork, info)
-            if (info /= 0) then
-            write(*,*) "Error in Eigenvalue decomposition of my info = ", info
-            stop
-            end if
-            deallocate(muy)
-            Yres(1:size(mx,1)) = vx(1:size(mx,1))
-            do i=1,size(mx,1)
-                  max_row=maxloc(abs(mx(:,i)),1)
-                  !yres(size(mx,1)+i)=(mx(max_row,i))**2
-                  !yres(size(mx,1)+i)=real(max_row)
-            enddo 
-
-            !Yres(size(mx,1)+1:2*size(mx,1)) = vy(1:size(my,1))
-      end subroutine
-
-      subroutine copy_2_upper(m)
-            implicit none
-            double precision, intent(inout) :: m(:,:)
-            integer :: i,j
-            ! copy the lower part of the matrix to the upper part
-            do i=1,size(m,1)
-                  do j=1,i-1
-                        m(j,i) = m(i,j)
-                  enddo
-            enddo 
-      end subroutine copy_2_upper
-
-      subroutine coppy_2_low(m)
-            implicit none
-            double precision, intent(inout) :: m(:,:)
-            integer :: i,j
-            ! copy the upper part of the matrix to the lower part
-            do i=1,size(m,1)
-                  do j=i+1,size(m,2)
-                        m(j,i) = m(i,j)
-                  enddo
-            enddo 
-      end subroutine coppy_2_low
-
-
-      !1 SUBROUTNE BLOCKS 
-      !! EACH BLOCK OF dIABTIC MATRIX 
-      SUBROUTINE block_diab(mx,my,Y,block)
-            implicit none
-            double precision, intent(inout):: Y(:)
-            double precision, intent(in) :: mx(ndiab,ndiab), my(ndiab,ndiab)
-            integer, intent(in) :: block
-            integer i,j,ii,nn
-            y=0.0d0 
-
-            select case (block)
-                  case(1)
-                        ! fill the first E1 block state 2 &3
-                        y(1)=mx(2,2)
-                        y(2)=mx(2,3)
-                        !y(3)=mx(3,2)
-                        y(4)=mx(3,3)
-                        !y(5)=my(2,2)
-                        !y(6)=my(2,3)
-                        !y(7)=my(3,2)
-                        !y(8)=my(3,3)
-
-                  case(2)
-                        ! fill the second E2 block state 4 & 5
-                        y(1)=mx(4,4)
-                        y(2)=mx(4,5)
-                        !y(3)=mx(5,4)
-                        y(4)=mx(5,5)
-                        y(5)=my(4,4)
-                        y(6)=my(4,5)
-                        !y(7)=my(5,4)
-                        y(8)=my(5,5)
-                  case(3)
-                        ! Filling the pseudo block E1 and E2
-                        y(1)=mx(2,4)
-                        y(2)=mx(2,5)
-                        y(3)=mx(3,4)
-                        y(4)=mx(3,5)
-                        y(5)=my(2,4)
-                        y(6)=my(2,5)
-                        y(7)=my(3,4)
-                        y(8)=my(3,5)
-                  case(4)
-                        ! filling the block of A2 coupling with E1
-                        y(1)=mx(1,2)
-                        y(2)=mx(1,3)
-                        y(3)=mx(2,1)
-                        y(4)=mx(3,1)
-                        !y(5)=my(1,2)
-                        !y(6)=my(1,3)
-                        !y(7)=my(2,1)
-                        !y(8)=my(3,1)
-                      
-                        
-                  case(5)
-                        ! couplinng A2 with E2
-                        Y(1)=mx(1,4)
-                        Y(2)=mx(1,5)
-                        !Y(3)=mx(4,1)
-                        !Y(4)=mx(5,1)
-                        Y(5)=my(1,4)
-                        Y(6)=my(1,5)
-                        !Y(7)=my(4,1)
-                        !Y(8)=my(5,1)
-                  case(6)
-                        ! Filling A only
-                        y(1)=mx(1,1)
-                        y(5)=my(1,1)
-
-                  case default
-                        write(*,*) "Error in block_diab subroutine, block not recognized"
-                        write(*,*) "The block is:", block
-                        stop
-            end select
-      end subroutine block_diab
-      subroutine ident(A)
-            implicit none
-            integer i,j
-            double precision,intent(inout)::A(:,:)
-            do i=1,size(A,1)
-            do j=1,size(A,1)
-            if (i==j) then 
-            A(i,j)=1.0d0
-            else
-            A(i,j)=0.0d0
-            endif
-            enddo 
-            enddo       
-      end subroutine
-
-
-      ! subroutine trasform the U matrix 
-      subroutine transform_U(U)
-            implicit none
-            double precision, intent(inout) :: U(ndiab,ndiab)
-          
-            double precision :: U_ref(ndiab,ndiab), V(ndiab,ndiab), E(nstat) 
-            integer i,max_row 
-            double precision:: dot_prod,q_ref(9)
-            logical,parameter:: dbg_sign =.true.
-
-            !q_ref= [1.000174,0.000000,0.000000,-0.503595,-0.872253,0.000000,-0.530624,0.919068,0.000000]
-            !call eval_surface(E,V,U_ref,q_ref,p) ! get the reference transformation matrix
-            do i=1,ndiab 
-                  max_row = maxloc(abs(U(:,i)),1)
-                  if (U(max_row,i) .lt. 0) then
-                        U(:,i) = -1*U(:,i)
-                  endif      
-            enddo
-            !dot_prod=dot_product(U(2:3,4),U_ref(2:3,4))
-            !if (dot_prod .lt. 0.0d0) then
-            !      U(:,4) = -1.0d0*U(:,4)
-            !endif
-           
-
-      end subroutine transform_U
-
-      subroutine write_type_calc(p,id_write)
-            ! Subroutine to write the type of calculation
-            implicit none
-            double precision, intent(in) :: p(:)
-            integer, intent(in) :: id_write
-            integer :: type_calc, blk
-            type_calc = int(p(pst(1,33)))
-            blk = int(p(pst(1,33)+1))
-
-            if (type_calc ==1) then 
-                  write(id_write,*) "Type of calculation: TRACE"
-            else if (type_calc ==2) then 
-                  write(id_write,*) "Type of calculation: EIGENVALUE"
-            else if (type_calc ==3) then 
-                  IF (blk == 1) then 
-                        write(id_write,*) "Type of calculation: E1 BLOCK"
-                  ELSE IF (BLK ==2) THEN 
-                        write(id_write,*) "Type of calculation: E2 BLOCK"
-                  ELSE IF (BLK ==3) THEN 
-                        write(id_write,*) "Type of calculation: Pseudo E1 and E2 BLOCK"
-                  ELSE IF (BLK ==4) THEN 
-                        write(id_write,*) "Type of calculation: COUPLING A2 with E1 BLOCK"
-                  ELSE IF (BLK ==5) THEN
-                        write(id_write,*) "Type of calculation: COUPLING A2 with E2 BLOCK"
-                  ELSE IF (BLK ==6) THEN 
-                        write(id_write,*) "Type of calculation: A2 ONLY"
-                  ELSE 
-                        write(id_write,*) "Type of calculation: Diabatic transformation with unknown block size", blk
-                  END IF
-      
-            else if (type_calc ==4) then 
-                  write(id_write,*) "Type of calculation: Full Diabatic Matrix"
-            else if (type_calc ==5) then 
-                  write(id_write,*) "Type of calculation: Transformation matrix U"
-            else 
-                  write(id_write,*) "Error in type of calculation:", type_calc
-                  stop
-            end if
-      END SUBROUTINE write_type_calc
-
-      !! subroutine for writting the transformtion matrix U
-      subroutine Transformation_mat(temp,v,y)
-            implicit none 
-            double precision, intent(in) :: temp(ndiab,ndiab), v(:)
-            double precision, intent(inout) :: y(:)
-            double precision :: U(ndiab,ndiab  )
-            integer i,j,ii
-            U(1:ndiab,1:ndiab) = temp(1:ndiab,1:ndiab)
-            !call transform_U(U,P)
-          
-            y=0.0d0
-            !y(1:4) = v(1:4) ! copy the first 4 elements of v to y
-            ii=1
-            do i=1,ndiab
-                  do j=1,ndiab
-                        y(ii) = U(i,j)
-                        ii=ii+1
-                  enddo
-            enddo
-            y(ii:30)=v(:)
-      end subroutine
-      ! compute the overlap between U matrix 
-      subroutine overlap(U_ref,U)
-            implicit none
-            double precision, intent(in):: U_ref(ndiab,ndiab)
-            double precision, intent(inout):: U(ndiab,ndiab)
-            double precision:: over
-            integer i
-            do i=1,ndiab
-
-                  over=dot_product(U_ref(:,i),U(:,i))
-                  if (over .lt. 0.0d0 ) then
-                        U(:,i)=-U(:,i)
-                  endif
-            enddo
-
-      end subroutine
-
-      
-end module 
diff --git a/src/model/model.f90 b/src/model/model.f90
deleted file mode 100644
index c970e80..0000000
--- a/src/model/model.f90
+++ /dev/null
@@ -1,507 +0,0 @@
-! Author: jnshuti
-! Created: 2025-10-03 14:09:49
-! Last modified: 2025-10-03 14:10:10 jnshuti
-! model for L-matrix of NO3 radical 
-
-module diab_mod
-  use accuracy_constants, only: dp, idp
-  use dim_parameter, only: ndiab, nstat, ntot,qn,pst
-  implicit none
-  private
-  public :: diab 
-  
-  contains 
-    subroutine diab(lx,ly,lz,n,x1,x2,p)
-        implicit none 
-        real(dp), intent(out),dimension(ndiab,ndiab):: lx,ly,lz 
-        real(dp), intent(in), dimension(qn):: x1,x2 
-        real(dp), intent(in),dimension(:):: p 
-        integer(idp),intent(in):: n
-
-        call Lx_diab(lx,x1,x2,p)
-        call Ly_diab(ly,x1,x2,p)
-        call Lz_diab(lz,x1,x2,p)
-
-    end subroutine diab 
-
-
-    subroutine Lx_diab(E,q,t,p)
-        implicit none 
-        real(dp),dimension(ndiab,ndiab), intent(out):: E 
-        real(dp),dimension(:),intent(in):: q,t
-        real(dp),dimension(:),intent(in):: p 
-        real(dp):: xs,ys,xb,yb,a,b 
-        real(dp):: v3_vec(8), v2(6)
-        integer(idp):: i,j,id
-
-
-
-        ! check the dimension of the matrix 
-        if (size(E,1) .ne. ndiab) then 
-            write(*,*) " Error in Lx_diab: wrong dimension of L matrix ", size(E,1)
-            stop
-        endif
-
-        ! rewrite the coordinate array q into symmetry adapted coordinates
-        call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
-        v2(1)=xs**2-ys**2
-        v2(2)=xb**2-yb**2
-        v2(3)=xs*xb-ys*yb
-        v2(4)=2*xs*ys 
-        v2(5)=2*xb*yb 
-        v2(6)=xs*yb+xb*ys
-
-        e = 0.0_dp 
-        id = 1
-        e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb ! 2 param
-        id=id+1   ! 2
-        e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb  ! 2 p
-        e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
-        id =id+1 ! 3
-        e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb   ! 2 p
-        e(5,5)=e(5,5)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
-
-        
-        id=id+1 ! 4
-        ! order 2
-
-        e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & ! 3 p
-                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
-        id =id+1 ! 5
-        e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
-                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
-        e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
-                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
-        id =id+1  ! 6
-        e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
-                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
-        e(5,5)=e(5,5)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
-                    +p(pst(1,id)+2)*(xs*xb-ys*yb)
-        ! W and Z term of E1
-        ! order 0
-        id=id+1  ! 7
-        e(2,2)=e(2,2)+p(pst(1,id))
-        e(3,3)=e(3,3)-p(pst(1,id))
-        !e(2,3)=e(2,3)
-
-        ! order 1
-        id=id+1  ! 8 ! 2 param
-        e(2,2)=e(2,2)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb
-        e(3,3)=e(3,3)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
-        e(2,3)=e(2,3)- p(pst(1,id))*ys -p(pst(1,id)+1)*yb 
-        ! order 2
-        id=id+1 ! 9  ! 3p 
-        do i=1,3
-        e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i)
-        e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i)
-        e(2,3)=e(2,3)+ p(pst(1,id)+(i-1))*v2(i+3)
-        enddo 
-        ! order 3
-        
-
-        ! try the testing of higher order terms 
-        !e(2,3)=e(2,3)- p(pst(1,id))*ys*ss +p(pst(1,id)+1)*ss*2*xs*ys 
-        
-
-
-
-        ! W and Z for E2
-        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        id=id+1  ! 10
-        e(4,4)=e(4,4)+p(pst(1,id))
-        e(5,5)=e(5,5)-p(pst(1,id))
-        e(4,5)=e(4,5)
-
-        ! order 1
-        id=id+1  !112 param 15
-        e(4,4)=e(4,4)+ p(pst(1,id))*xs+p(pst(1,id)+1)*xb  
-        e(5,5)=e(5,5)- (p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
-        e(4,5)=e(4,5)- p(pst(1,id))*ys-p(pst(1,id)+1)*yb 
-        ! order 2
-        id=id+1 ! 12 ! 3p 
-        do i=1,3
-        e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i)
-        e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i)
-        e(4,5)=e(4,5)+ p(pst(1,id)+(i-1))*v2(i+3)
-        enddo 
-
-        ! make the dipole E = b* E
-
-        e = b * e 
-
-        ! E1 X E2 
-        ! WW and ZZ 
-        
-
-        id =id+1 ! 13
-        e(2,4)=e(2,4)+p(pst(1,id))*b
-        e(3,5)=e(3,5)-p(pst(1,id))*b
-
-        ! ORDER 1
-        id=id+1 ! 14  ! 6 parama
-        e(2,4)=e(2,4)+b*((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*xb)
-        e(3,5)=e(3,5)+b*((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
-        e(2,5)=e(2,5)+b*((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
-        e(3,4)=e(3,4)+b*((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
-        ! order 2
-        id=id+1   ! 15
-
-        do i=1,3 !  param
-        e(2,4)=e(2,4)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
-        e(3,5)=e(3,5)+b*(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
-        e(2,5)=e(2,5)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
-        e(3,4)=e(3,4)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i+3)
-        enddo 
-        
-        ! pseudo A2 & E1
-        ! ##################################################
-        !###################################################
-        ! order 0
-        id=id+1 ! 1 param ! 16
-
-        e(1,3)=e(1,3)+b*(p(pst(1,id)))
-        
-        ! order 1 
-        id = id +1 ! 17 
-        e(1,2)=e(1,2)-b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
-        e(1,3)=e(1,3)+b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
-        
-        ! order 2 
-        id=id+1 ! 18
-
-        e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
-        e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) &
-                        +p(pst(1,id)+2)*(xs*xb-ys*yb))
-
-        ! COUPLING OF A2 WITH E2
-
-        !##########################################################################################################
-                    
-        ! order 0
-        id =id+1 !19
-        e(1,5)=e(1,5)+p(pst(1,id))
-        
-        ! order 1 
-        id = id +1 ! 20 
-        e(1,4)=e(1,4)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
-        e(1,5)=e(1,5)+(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
-        
-        ! order 2 
-        id=id+1 ! 21
-
-        e(1,4)=e(1,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys)
-        e(1,5)=e(1,5)+p(pst(1,id))*(xs**2-ys**2) + p(pst(1,id)+1)*(xb**2-yb**2) &
-                        +p(pst(1,id)+2)*(xs*xb-ys*yb)
-
-        e(1,4:5) = b* e(1,4:5)
-        
-
-        call copy_2_lower_triangle(e)
-    end subroutine Lx_diab
-
-        ! Ly matrix
-    subroutine Ly_diab(e,q,t,p)
-        implicit none 
-        real(dp),dimension(ndiab,ndiab), intent(out):: e
-        real(dp),dimension(:),intent(in):: q,t
-        real(dp),dimension(:),intent(in):: p 
-        real(dp):: xs,ys,xb,yb,a,b 
-        real(dp):: v2(6)
-        integer(idp):: i,j,id
-
-        ! check the dimension of the matrix 
-        if (size(e,1) .ne. ndiab) then 
-            write(*,*) " Error in Ly_diab: wrong dimension of L matrix ", size(e,1)
-            stop
-        endif
-
-        ! rewrite the coordinate array q into symmetry adapted coordinates
-        call rewrite_coord(q,a,xs,ys,xb,yb,b,1) 
-        v2(1)=xs**2-ys**2
-        v2(2)=xb**2-yb**2
-        v2(3)=xs*xb-ys*yb
-        v2(4)=2*xs*ys 
-        v2(5)=2*xb*yb 
-        v2(6)=xs*yb+xb*ys
-
-        e = 0.0_dp 
-        
-        ! V-term 
-
-        id=1  ! 1
-        ! order 1
-        e(1,1)=e(1,1)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
-        id=id+1  ! 2
-        e(2,2)=e(2,2)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
-        e(3,3)=e(3,3)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
-        id =id+1  ! 3
-        e(4,4)=e(4,4)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
-        e(5,5)=e(5,5)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb
-
-        id=id+1  ! 4b*(
-        e(1,1)=e(1,1)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
-        id =id+1 ! 5
-        e(2,2)=e(2,2)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
-        e(3,3)=e(3,3)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
-        id=id+1 ! 6
-        e(4,4)=e(4,4)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
-        e(5,5)=e(5,5)-(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+p(pst(1,id)+2)*(xs*yb+xb*ys))
-         
-        
-        !  W and Z of E1
-        ! order 0
-        id=id+1 ! 7
-        e(2,3)=e(2,3)+p(pst(1,id))
-        ! order 1
-        id=id+1 ! 8
-        e(2,2)=e(2,2)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb 
-        e(3,3)=e(3,3)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb 
-        e(2,3)=e(2,3)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
-        ! order 2
-        id=id+1 ! 9
-        do i=1,3
-        e(2,2)=e(2,2)+p(pst(1,id)+(i-1))*v2(i+3)
-        e(3,3)=e(3,3)-p(pst(1,id)+(i-1))*v2(i+3)
-        e(2,3)=e(2,3)-p(pst(1,id)+(i-1))*v2(i) 
-        enddo 
-        
-
-        !! W and Z of E2
-        !!!!!!!!!!!!!!!!!!!!!!!!!
-        !!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-        ! order 0
-        id=id+1 ! 10
-        e(4,5)=e(4,5)+p(pst(1,id))
-        ! order 1
-        id=id+1 ! 11
-        e(4,4)=e(4,4)-p(pst(1,id))*ys -p(pst(1,id)+1)*yb 
-        e(5,5)=e(5,5)+p(pst(1,id))*ys+ p(pst(1,id)+1)*yb 
-        e(4,5)=e(4,5)-p(pst(1,id))*xs -p(pst(1,id)+1)*xb
-        ! order 2
-        id=id+1 ! 12
-        do i=1,3
-        e(4,4)=e(4,4)+p(pst(1,id)+(i-1))*v2(i+3)
-        e(5,5)=e(5,5)-p(pst(1,id)+(i-1))*v2(i+3)
-        e(4,5)=e(4,5)-p(pst(1,id)+(i-1))*v2(i) 
-        enddo 
-
-        ! PSEUDO JAHN-TELLER E1 AND E2
-
-        e = b* e 
-
-        !ORDER 0 
-        id=id+1 ! 13
-
-        e(2,5)=e(2,5)+p(pst(1,id))
-        e(3,4)=e(3,4)+p(pst(1,id))
-        ! order 1
-        
-        id=id+1 ! 14
-        e(2,4)=e(2,4)+((p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*yb)
-        e(3,5)=e(3,5)+((p(pst(1,id))+p(pst(1,id)+1)+p(pst(1,id)+2))*ys+(p(pst(1,id)+3)+p(pst(1,id)+4)+p(pst(1,id)+5))*yb)
-        e(2,5)=e(2,5)+((-p(pst(1,id))+p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(-p(pst(1,id)+3)+p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
-        e(3,4)=e(3,4)+((p(pst(1,id))-p(pst(1,id)+1)-p(pst(1,id)+2))*xs+(+p(pst(1,id)+3)-p(pst(1,id)+4)-p(pst(1,id)+5))*xb)
-
-        ! order 2
-        id=id+1 ! 15 
-
-        e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
-        e(3,5)=e(3,5)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i+3)
-        e(2,5)=e(2,5)+(-p(pst(1,id)+(i-1))+p(pst(1,id)+(i+2))-p(pst(1,id)+(i+5)))*v2(i)
-        e(3,4)=e(3,4)+(-p(pst(1,id)+(i-1))-p(pst(1,id)+(i+2))+p(pst(1,id)+(i+5)))*v2(i)
-
-        ! no order 3
-        !!!!!!!!!!!!!!!!
-
-        ! the coupling A2 & E1
-        ! #####################
-        ! order 0
-
-        id=id+1 ! 16 
-        e(1,2)=e(1,2)+(p(pst(1,id)))
-        ! order 1
-
-        id=id+1 ! 17 
-        e(1,2)=e(1,2)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
-        e(1,3)=e(1,3)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
-
-        ! order 2
-        id=id+1 !18 
-        e(1,2)=e(1,2)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
-                        +p(pst(1,id)+2)*(xs*xb-ys*yb))
-        e(1,3)=e(1,3)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ &
-                        +p(pst(1,id)+2)*(xs*yb+xb*ys))
-
-        ! COUPLING OF A2 WITH E2
-        !#######################################################################################
-        !###############################################################################
-        ! order 0 
-
-        id = id+1 !19
-        e(1,4)=e(1,4)+p(pst(1,id))
-        ! order 1
-
-        id=id+1 ! 20
-        e(1,4)=e(1,4)-(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
-        e(1,5)=e(1,5)-(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
-
-        ! order 2
-        id=id+1 ! 21 
-        e(1,4)=e(1,4)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
-                        +p(pst(1,id)+2)*(xs*xb-ys*yb))
-        e(1,5)=e(1,5)+(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)+ &
-                        p(pst(1,id)+2)*(xs*yb+xb*ys))
-
-        !write(*,*)'idy=',id
-        e(1:4,5) = b * e(1:4,5)
-
-
-        call copy_2_lower_triangle(e)
-
-    end subroutine Ly_diab
-         ! Lz matrix
-    subroutine Lz_diab(e,q,t,p)
-        implicit none 
-        real(dp),dimension(ndiab,ndiab), intent(out):: e 
-        real(dp),dimension(:),intent(in):: q,t
-        real(dp),dimension(:),intent(in):: p 
-        real(dp):: xs,ys,xb,yb,a,b 
-        real(dp):: v2(6)
-        integer(idp):: i,j,id 
-
-        ! check the dimension of the matrix 
-        if (size(e,1) .ne. ndiab) then 
-            write(*,*) " Error in Lz_diab: wrong dimension of e matrix ", size(e,1)
-            stop
-        endif
-        call rewrite_coord(q,a,xs,ys,xb,yb,b,1)
-
-
-        e = 0.0_dp
-        ! id for lz 
-        id =  22 ! has to be 
-        ! the diagonal terms 
-
-        ! the v-term is 0th order and 3rd order.
-        ! There is no zeroth order for diagonal
-
-        ! w and z of E''
-        ! order 1 
-        id = id  ! 22
-        e(2,2) = e(2,2) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb  
-        e(3,3) = e(3,3) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
-        e(2,3) = e(2,3) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb
-
-        ! order 2
-        id = id +1 ! 23
-        do i =1,3 
-            e(2,2) = e(2,2) + p(pst(1,id)+(i-1))*v2(i+3) 
-            e(3,3) = e(3,3) - p(pst(1,id)+(i-1))*v2(i+3)
-            e(2,3) = e(2,3) + p(pst(1,id)+(i-1))*v2(i)
-        enddo 
-
-        ! W and Z of E'
-        ! order 1 
-        
-        id = id +1 ! 24
-        e(4,4) = e(4,4) + p(pst(1,id))*ys + p(pst(1,id)+1)*yb  
-        e(5,5) = e(5,5) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb
-        e(4,5) = e(4,5) - p(pst(1,id))*xs -p(pst(1,id)+1)*xb
-
-        ! order 2
-        id = id +1 ! 25
-         do i =1,3 
-            e(4,4) = e(4,4) + p(pst(1,id)+(i-1))*v2(i+3) 
-            e(5,5) = e(5,5) - p(pst(1,id)+(i-1))*v2(i+3)
-            e(4,5) = e(4,5) + p(pst(1,id)+(i-1))*v2(i)
-        enddo 
-        ! the coupling 
-        ! Pseudo of E' and E'' 
-        ! it must have odd power of b
-
-        id = id +1 !26
-        ! order 0 
-        e(2,4) = e(2,4) 
-        e(3,5) = e(3,5) 
-        e(2,5) = e(2,5) + b*(p(pst(1,id)))
-        e(3,4) = e(3,4) - b*(p(pst(1,id)))
-
-        ! order 1 
-        id = id +1 !27
-        e(2,4) = e(2,4) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) 
-        e(3,5) = e(3,5) + b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
-        e(2,5) = e(2,5) - b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
-        e(3,4) = e(3,4) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*xb)
-        ! order 2 
-        id = id +1 !28
-        do i=1,3 
-            e(2,4) = e(2,4) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
-            e(3,5) = e(3,5) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
-            e(2,5) = e(2,5) + b*(p(pst(1,id)+(i-1)))*v2(i)
-            e(3,4) = e(3,4) - b*(p(pst(1,id)+(i-1)))*v2(i)
-        enddo 
-        ! no third order 
-
-        ! the coupling between A2' and E''
-        ! order 1 
-        id = id +1 !29
-        e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)*xb))
-        e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)*yb))
-
-
-        id = id +1 !30
-        ! order 2 
-        do i=1,3
-            e(1,2) = e(1,2) + b*(p(pst(1,id)+(i-1)))*v2(i)
-            e(1,3) = e(1,3) + b*(p(pst(1,id)+(i-1)))*v2(i+3)
-        enddo 
-
-        ! the coupling of A2' and E'
-
-        ! order 1 
-        id = id +1 !31
-        e(1,2) = e(1,2) + (p(pst(1,id))*xs + p(pst(1,id)*xb))
-        e(1,3) = e(1,3) - (p(pst(1,id))*ys + p(pst(1,id)*yb))
-        
-        id = id +1 ! 32
-        ! order 2 
-        do i=1,3
-            e(1,2) = e(1,2) + (p(pst(1,id)+(i-1)))*v2(i)
-            e(1,3) = e(1,3) + (p(pst(1,id)+(i-1)))*v2(i+3)
-        enddo 
-        call copy_2_lower_triangle(e)
-
- 
-
-    end subroutine Lz_diab
-
-    subroutine rewrite_coord(q,a,xs,ys,xb,yb,b,start)
-        implicit none 
-        real(dp),dimension(:),intent(in):: q
-        real(dp),intent(out):: xs,ys,xb,yb,a,b
-        integer(idp),intent(in):: start
-        integer(idp):: i,j
-
-        a= q(start)
-        xs = q(start+1)
-        ys = q(start+2)
-        xb = q(start+3)
-        yb = q(start+4)
-        b = q(start+5)
-    end subroutine rewrite_coord
-
-    subroutine copy_2_lower_triangle(mat)
-        real(dp), intent(inout) :: mat(:, :)
-        integer :: m, n
-    !     write lower triangle of matrix symmetrical
-        do n=1,size(mat,1)
-            do m=n,size(mat,1)
-                mat(m,n)=mat(n,m)
-            enddo 
-        enddo
-   end subroutine copy_2_lower_triangle
-
-
-end module diab_mod
diff --git a/src/model/nnparams.incl b/src/model/nnparams.incl
deleted file mode 100644
index 635c33c..0000000
--- a/src/model/nnparams.incl
+++ /dev/null
@@ -1,43 +0,0 @@
-!**** 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)
diff --git a/src/model/surface_mod.f90 b/src/model/surface_mod.f90
deleted file mode 100644
index b9b7315..0000000
--- a/src/model/surface_mod.f90
+++ /dev/null
@@ -1,85 +0,0 @@
-module surface_mod
-    use accuracy_constants, only: dp
-    implicit none
-    private
-    public eval_surface
-    contains
- subroutine eval_surface(e, w, u, x1)
-   
-    use accuracy_constants, only: dp, idp
-    use dim_parameter, only: ndiab
-    implicit none
-    real(dp), dimension(:, :), intent(out) :: w, u
-    real(dp), dimension(:), intent(out) :: e
-    real(dp), dimension(:), intent(in) :: x1
-    real(dp), allocatable, dimension(:, :) :: Mat
-
-    ! debug parameter
-
-    logical, parameter:: dbg=.false.
-    integer(kind=idp):: i,j
-      !     lapack variables
-      integer(kind=idp), parameter :: lwork = 1000
-      real(kind=dp) work(lwork)
-      integer(kind=idp) info
-      
-
-      !write(*,*)"# Calling the potential routine "
-      call init_pot_para
-      call potentialno35s(W,X1)
-      
-      allocate (Mat, source=w)
-      call dsyev('V', 'U', ndiab, Mat, ndiab, e, work, lwork, info)
-      if( info .ne. 0) then 
-         write(*,*) " Error in eigenvalues decomposition routine of potential info=", info 
-         stop 
-      endif 
-      u(:, :) = Mat(:, :)
-      deallocate (Mat)
-      
-      if (dbg) then 
-         do i=1,ndiab 
-            
-            write(19,99) e(i),(U(i,j),j=1,ndiab)
-         enddo 
-         write(19,*)""
-      endif 
- 99 format(2x,f16.8,2X,5f16.8)
- 
- end subroutine eval_surface
- 
-! subroutine init_surface(p)
-!    use dim_parameter, only: ndiab, nstat, ntot, nci ,qn
-!    use parameterkeys, only: parameterkey_read
-!    use fileread_mod, only: get_datfile, internalize_datfile
-!    use io_parameters, only: llen
-!    use accuracy_constants, only: dp
-!    implicit none
-!    real(dp), dimension(:), allocatable, intent(out) :: p
-!    character(len=llen), allocatable, dimension(:) :: infile
-! 
-!    qn = 9
-!    ndiab = 4
-!    nstat = 4
-!    nci = 4
-!    ntot = ndiab + nstat + nci
-! 
-!    block
-!       character(len=:),allocatable :: datnam
-!       integer :: linenum
-!       !get parameter file
-!       call get_datfile(datnam)
-!       !internalize datfile
-!       call internalize_datfile(datnam, infile, linenum, llen)
-!    end block
-! 
-!    !read parameters from file
-!    block
-!       real(dp), dimension(:), allocatable :: p_spread
-!       integer,dimension(:),allocatable :: p_act
-!       integer :: npar
-!       real(dp), parameter :: facspread = 1.0_dp, gspread = 1.0_dp
-!       call parameterkey_read(infile, size(infile, 1), p, p_act, p_spread, npar, gspread, facspread)
-!    end block
-! end subroutine init_surface
- end module surface_mod
diff --git a/src/model/weight.f b/src/model/weight.f
deleted file mode 100644
index 99164c5..0000000
--- a/src/model/weight.f
+++ /dev/null
@@ -1,50 +0,0 @@
-!     <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
diff --git a/src/model/write.f b/src/model/write.f
deleted file mode 100644
index 700fd30..0000000
--- a/src/model/write.f
+++ /dev/null
@@ -1,763 +0,0 @@
-      module write_mod
-      implicit none
-!     unit conversion
-      double precision ,parameter :: h2icm = 219474.69d0
-      double precision, parameter :: au2Debye = 2.541746d0
-      character(len=250), parameter :: sep_line = '(250("-"))'
-      character(len=250), parameter :: block_line = '(250("="))'
-
-      contains
-
-!     <Subroutine for writing the Output
-      subroutine write_output
-     >     (q,x1,x2,y,wt,par,p_act,p_spread,nset,npar,
-     >     flag,lauf)
-      use adia_mod, only: adia
-      use dim_parameter,only: qn,ntot,numdatpt,ndiab
-      use ctrans_mod,only: ctrans
-      implicit none
-!     IN: variables
-      integer lauf
-      integer flag              !< 0= initial output 1=fit not converged 2= Fit Converged, 3= max iteration reached
-      integer npar,nset
-      double precision par(npar,nset),p_spread(npar)
-      integer p_act(npar)
-      double precision q(qn,numdatpt),x1(qn,numdatpt),x2(qn,numdatpt)
-      double precision y(ntot,numdatpt),wt(ntot,numdatpt)
-
-!     INTERNAL: Variables
-      integer,parameter  :: id_out = 20 , std_out = 6
-      integer pt
-      integer i, id_print
-      double precision, allocatable ::  ymod(:,:)
-      double precision, allocatable ::  ew(:,:)
-      double precision, allocatable ::  ev(:,:,:)
-
-      logical skip
-
-      allocate(ymod(ntot,numdatpt))
-      allocate(ew(ndiab,numdatpt))
-      allocate(ev(ndiab,ndiab,numdatpt))
-
-      skip=.false.
-
-!     get Model Outputs for all geometries for current best parameter set par(:,1)
-      do pt=1,numdatpt
-         call adia(pt,par(1:npar,1),npar,ymod(1:ntot,pt),
-     >        ew(1:ndiab,pt),ev(1:ndiab,1:ndiab,pt),skip)
-         call ctrans(q(:,pt),x1(:,pt),x2(:,pt))
-      enddo
-
-!     Initial write print everything you want to see before the fit and return
-      if(flag.eq.0) then
-         call print_parameterstate(std_out,par(:,1),p_act,npar)
-         call print_ErrorSummary(std_out,y,ymod,wt)
-!        print Data into the plotfiles
-         return
-      endif
-!     open output files for individual makro iterations
-      call open_outfile(id_out,lauf)
-!     print Data into the plotfiles
-      call print_plotfiles(x1,y,wt,ymod)
-
-!     print Genetic output into files
-      do i=1, 2
-         if (i.eq.1) then
-            id_print= std_out
-         else
-            id_print= id_out
-         endif
-         write(id_print,'("Writing Iteration: ",i4)') lauf
-         write(id_print,block_line)
-!        write data information only in outfile
-         if(i.eq.2) then
-            call print_data(id_print,x1,y,ymod,wt)
-            call print_Set_Errors(id_print,y,ymod,wt)
-         endif
-         call print_parameterblock
-     >        (id_print,par(:,1),p_act,p_spread,npar)
-         call print_ErrorSummary(id_print,y,ymod,wt)
-
-      enddo
-
-      call print_fortranfile(par(:,1),npar)
-      
-      ! write the type of calc at the end of the output
-
-
-      close (id_out)
-      deallocate(ymod,ev,ew)
-      end subroutine
-!----------------------------------------------------------------------------------------------------
-!     <subroutine for scan seperated Error analysis>
-      subroutine print_Set_Errors(id_out,y, ymod, wt)
-      use io_parameters,only: llen
-      use dim_parameter,only: ndata,nstat,ntot,numdatpt,sets
-      integer , intent(in) :: id_out
-      double precision, intent(in) ::  y(ntot,numdatpt),
-     >     ymod(ntot,numdatpt), wt(ntot,numdatpt)
-      integer :: set, setpoint, pt
-      double precision :: Set_rms(sets,ntot), Set_num(sets,ntot)
-      double precision :: Total_rms, Total_Energy_rms,Energy_rms(nstat)
-      character(len=llen) fmt
-      write(id_out,'(A)') 'Errors in icm for individual Sets' //
-     >     '(specified by sets: and npoints:)'
-      write(id_out,'(A5,3A16)')'Set','Total',
-     >     'Total_Energy', 'Energy[nstat]'
-      write(id_out,sep_line)
-      write(fmt,'("(I5,2f16.1,",I2,"f16.1)")') nstat
-      Set_rms = 0.d0
-      pt = 0
-      do set=1, sets
-         do setpoint=1, ndata(set)
-            pt = pt + 1
-            where(wt(:,pt) > 0.d0)
-               Set_rms(set,:) = Set_rms(set,:)+(ymod(:,pt)-y(:,pt))**2
-               Set_num(set,:) = Set_num(set,:) + 1
-            end where
-         enddo
-         Total_rms
-     >        = dsqrt(sum(Set_rms(set,:))
-     >        / (sum(Set_num(set,:))))
-         Total_Energy_rms
-     >        = dsqrt(sum(Set_rms(set,1:nstat))
-     >        / (sum(Set_num(set,1:nstat))))
-         Energy_rms(1:nstat)
-     >        = dsqrt(Set_rms(set,1:nstat)
-     >        / (Set_num(set,1:nstat)))
-         write(id_out,fmt) set, Total_rms*h2icm, Total_Energy_rms*h2icm,
-     >        Energy_rms(1:nstat)*h2icm
-      enddo
-      write(id_out,block_line)
-      write(id_out,*) ''
-      end subroutine print_Set_Errors
-
-!----------------------------------------------------------------------------------------------------
-!     <subroutine for printing the parameter and the pst vector in fortran readable style for including the fitted parameters in other programs
-      subroutine print_fortranfile(p,npar)
-      use io_parameters,only: maxpar_keys
-      use dim_parameter,only: pst
-      implicit none
-!     IN: variables
-      integer npar
-      double precision p(npar)
-!     INTERNAL: variables
-      integer i
-      integer, parameter :: id_out = 49
-      character(len=32), parameter :: fname ='fit_genric_bend_no3.f90'
-
-      open(id_out,file=fname)
-
- 30   format(6x,A2,i3,A2,d18.9)
- 31   format(6x,A6,i3,A2,i3)
-      
-      write(id_out,'(2X,A)') "Module dip_param"
-      write(id_out,'(5X,A)') "IMPLICIT NONE"
-      write(id_out,'(5X,A,I0)') "Integer,parameter :: np=",npar 
-      write(id_out,'(5X,A,I0,A)') "Double precision :: p(",npar,")"
-      write(id_out,'(5X,A,I0,A)') "integer :: pst(2,",maxpar_keys,")"
-      write(id_out,'(5X,A)') "contains"
-      write(id_out,*)''
-  
-      write (id_out,'(5x,a)') "SUBROUTINE init_dip_planar_data()"
-      write (id_out,'(8X,A)') "implicit none"
-      do i=1,npar
-         write(id_out,30) 'p(',i,')=',p(i)
-      enddo
-      do i=1,maxpar_keys
-         write(id_out,31) 'pst(1,',i,')=',pst(1,i)
-         write(id_out,31) 'pst(2,',i,')=',pst(2,i)
-      enddo
-      
-
-      write(id_out,"(A)") "End SUBROUTINE init_dip_planar_data"
-      write(id_out,"(A)") "End Module dip_param"
-
-      close(id_out)
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <subroutine print_ErrorSummary: calculates the rms errros and prints them in the corresponding file
-      subroutine print_ErrorSummary(id_out,y,ymod,wt)
-      use dim_parameter,only: nstat,rms_thr,ntot,numdatpt
-      use io_parameters,only: llen
-      implicit none
-!     IN: variables
-      integer id_out
-      double precision y(ntot,numdatpt),ymod(ntot,numdatpt)
-      double precision wt(ntot,numdatpt)
-!     INTERNAL: variables
-!     Counter and RMS variables
-      double precision  Cut_thr(nstat)
-      double precision  Output_rms(ntot),Cut_rms(nstat),Weighted_rms
-      integer Output_num(ntot),Cut_num(nstat)
-      double precision  Weighted_wt
-      double precision  Total_rms,Total_Weighted_rms
-      double precision  Total_Energie_rms,Total_State_rms(nstat)
-      double precision  Cut_Energie_rms, Cut_State_rms(nstat)
-
-      ! Variables for computing the NRMSE 
-      !double precision:: ymean(ntot),ysum(ntot),NRMSE 
-
-!     loop control
-      integer j,pt
-
-!     Fabian
-      character(len=llen) fmt
-!     initialize RMS variables
-      Output_rms(1:ntot) = 0.d0
-      Output_num(1:ntot) = 0
-      Weighted_rms  = 0.d0
-      Weighted_wt = 0.d0
-      Cut_rms(1:nstat)= 0.d0
-      Cut_num(1:nstat)= 0
-
-!     Define Threshold for Cut_* RMS Values
-      Cut_thr(1:nstat) = rms_thr(1:nstat)
-!     SUMM!
-!     Loop over all Datapoints
-      do pt=1,numdatpt
-!        get unweighted rms for each output value and count their number
-         do j=1,ntot
-            if(wt(j,pt).gt.0.d0) then
-               Output_rms(j) = Output_rms(j) +
-     >              (ymod(j,pt)-y(j,pt))**2
-               Output_num(j)=Output_num(j) + 1
-            endif
-         enddo
-!        get the unweighted rms under the given threshold and count their number
-         do j=1,nstat
-            if(wt(j,pt).gt.0.d0) then
-               if(y(j,pt).le.Cut_thr(j)) then
-                  Cut_rms(j) = Cut_rms(j) +
-     >                 (ymod(j,pt)-y(j,pt))**2
-                  Cut_num(j) = Cut_num(j) + 1
-               endif
-            endif
-         enddo
-!        get the weighted rms over all output values
-         Weighted_rms = Weighted_rms +
-     >        sum(((ymod(1:ntot,pt)-y(1:ntot,pt))**2)
-     >        *(wt(1:ntot,pt)**2))
-         Weighted_wt = Weighted_wt + sum(wt(1:ntot,pt)**2)
-      enddo
-
-!     NORM!
-!     TOTAL RMS:
-!     unweighted
-      Total_rms =
-     >     dsqrt(sum(Output_rms(1:ntot)) /(sum(Output_num(1:ntot))))
-
-!     Weighted
-      Total_Weighted_rms = dsqrt(Weighted_rms/Weighted_wt)
-
-!     unweighted, considering only first nstat values
-      Total_Energie_rms =
-     >     dsqrt(sum(Output_rms(1:nstat)) /(sum(Output_num(1:nstat))))
-
-!     unweighted,for each of the first nstat values separatly
-      Total_State_rms(1:nstat) =
-     >     dsqrt(Output_rms(1:nstat) / Output_num(1:nstat))
-
-!     unweighted,first nstat values only counting points under given threshold
-      Cut_Energie_rms =
-     >     dsqrt(sum(Cut_rms(1:nstat)) /(sum(Cut_num(1:nstat))))
-
-!     unweighted,each nstat values seperatly only counting points under threshold
-      Cut_State_rms(1:nstat) =
-     >     dsqrt(Cut_rms(1:nstat)/Cut_num(1:nstat))
-
-!     WRITE!
-!     make the actual writing into the file
-      write(id_out,39)
-      write(id_out,40)
-      write(id_out,41) Total_rms, Total_rms*au2Debye!Total_rms*h2icm
-      write(id_out,42) sum(Output_num(1:ntot))
-      write(id_out,43) Total_Weighted_rms, Total_Weighted_rms*h2icm
-      write(id_out,44) Weighted_wt
-      write(id_out,45) Total_Energie_rms, Total_Energie_rms*h2icm
-      write(id_out,42) sum(Output_num(1:nstat))
-      write(fmt,'("(A,10x,A,",I2,"f8.1)")') nstat
-      write(id_out,fmt) '#','State resolved RMS(icm): ',
-     $     Total_State_rms(1:nstat)*h2icm
-      write(fmt,'("(A,10x,A,",I2,"i8)")') nstat
-      write(id_out,fmt) '#','No. of Points per State: ',
-     $     Output_num(1:nstat)
-      write(id_out,51)
-
-!     write the errors under a given threshold if there were any points
-      if(any(Cut_num(1:nstat).gt.0)) then
-         write(id_out,48) Cut_Energie_rms, Cut_Energie_rms*h2icm
-         write(id_out,42) sum(Cut_num(1:nstat))
-
-         write(fmt,'("(A,10x,A,",I2,"f8.1,A)")') nstat
-         write(id_out,fmt) '#','Red. State resolved RMS: ',
-     $        Cut_State_rms(1:nstat)*h2icm,' icm'
-         write(fmt,'("(A,10x,A,",I2,"i8)")') nstat
-         write(id_out,fmt) '#','No. of Points per State: ',
-     $        Cut_num(1:nstat)
-         write(fmt,'("(A,10x,A,",I2,"f8.1,A)")') nstat
-         write(id_out,fmt) '#','Threshold per State:     ',
-     $        Cut_thr(1:nstat)*h2icm,' icm above Reference Point.'         
-
-      endif
-      write(id_out,39)
-
-!     FORMAT! specifications for the writing
- 39   format(250('#'))
- 40   format('#',10x,'ERROR SUMMARY: ')
- 41   format('#',10x,'Total RMS:            ',g16.8, '(',g16.8,
-     >     ' Debye)')
- 42   format('#',10x,'No. of Points:        ',i10)
- 43   format('#',10x,'Total weighted RMS:   ',g16.8, '(',f8.1,' icm)')
- 44   format('#',10x,'Sum of point weights: ',f16.8)
- 45   format('#',10x,'Total Energie RMS:    ',g16.8, '(',f8.1,' icm)')
-
- 48   format('#',10x,'Red. Energie RMS:        ',g16.8,'(',f8.1,' icm)')
- 51   format('#')
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-      subroutine print_plotfiles(x,y,wt,ymod)
-      use dim_parameter,only: ndata,sets,qn,ntot,numdatpt,plot_coord
-      implicit none
-!     IN: variables
-      double precision x(qn,numdatpt),y(ntot,numdatpt)
-      double precision wt(ntot,numdatpt), ymod(ntot,numdatpt)
-!     INTERNAL: variables
-      integer sstart,ssend,set,id_plot
-
-!     Initialize position pointer
-      ssend=0
-!     loop over datasets and print the plotfiles
-      do set=1 ,sets
-         if(ndata(set).eq.0) cycle
-         id_plot=50+set
-         call open_plotfile(id_plot,set)
-         write(id_plot,'(A)') '# -*- truncate-lines: t -*-'
-!        get start and end point of each set
-         sstart=ssend+1
-         ssend=ssend+ndata(set)
-         if (plot_coord(set).eq.0) then
-            call print_plotwalk(x(:,sstart:ssend),y(:,sstart:ssend),
-     >        wt(:,sstart:ssend),ymod(:,sstart:ssend),
-     >           ndata(set),id_plot,set)
-         else
-            call print_plotcoord(plot_coord(set),
-     >           x(:,sstart:ssend),y(:,sstart:ssend),
-     >           wt(:,sstart:ssend),ymod(:,sstart:ssend),
-     >           ndata(set),id_plot,set)
-         endif
-         close(id_plot)
-      enddo
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-      subroutine print_plotwalk(x,y,wt,ymod,npt,id_plot,set)
-      use dim_parameter,only: qn,ntot
-      use io_parameters,only: llen
-      implicit none
-!     IN: variables
-      integer id_plot,npt,set
-      double precision x(qn,npt),y(ntot,npt),ymod(ntot,npt),wt(ntot,npt)
-!     INTERNAL: variables
-      double precision xdiff(qn),walktime
-      double precision walknorm
-!     loop control
-      integer i,j
-
-      character(len=llen) fmt
-      j=ntot-1
-
-      call print_plotheader(id_plot,0,npt,set)
-
-      call getwalknorm(x,walknorm,npt)
-      walktime = 0.d0
-      do i=1,npt
-         if(i.gt.1) then
-            xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1)
-            walktime = walktime + dsqrt(sum(xdiff(1:qn)**2))/walknorm
-         endif
-         write(id_plot,"(ES16.8,*(3(ES16.8),:))")
-     >        walktime ,ymod(:,i),y(:,i),(wt(:,i))
-      enddo
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-      subroutine print_plotcoord(coord,x,y,wt,ymod,npt,id_plot,set)
-      use dim_parameter,only: qn,ntot
-      use io_parameters,only: llen
-      implicit none
-!     IN: variables
-      integer, intent(in) :: id_plot,npt,set,coord
-      double precision, intent(in) :: x(qn,npt),y(ntot,npt)
-      double precision, intent(in) :: ymod(ntot,npt),wt(ntot,npt)
-!     loop control
-      integer i
-
-      call print_plotheader(id_plot,coord,npt,set)
-      do i=1,npt
-!         write(id_plot,"(ES16.8,*(3(ES16.8),:))")
-!     >        x(coord,i), ymod(:,i),y(:,i),(wt(:,i))
-         write(id_plot,"(2ES16.8,*(3(ES16.8),:))")
-     >        x(coord,i), x(coord+1,i),y(:,i)
-      enddo
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-      subroutine print_plotheader(id_plot,coord,npt,set)
-      use dim_parameter,only: qn,ntot
-      use io_parameters,only: llen
-      implicit none
-      integer, intent(in) :: id_plot,npt,set,coord
-
-      character(len=llen) fmt
-
-      write(id_plot,'("#SET: ",i5)') set
-      write(id_plot,'("#OUTPUT VALUES",i4)') ntot
-      write(id_plot,'("#DATA POINTS: ",i4)') npt
-      if (coord.le.0) then
-      write(id_plot,'("#t(x) = WALK")')
-      else
-         write(id_plot,'("#t(x) = x(",I0,")")') coord
-      endif
-      write(id_plot,'("#UNIT: hartree")')
-      write(id_plot,'()')
-      write(id_plot,'("#",A15)',advance='no') "t(x)"
-      write(fmt,'("(3(7X,A9,",I3,"(16x)))")') ntot-1
-      write(id_plot,fmt) 'ymod(p,x)','y(x)  ','wt(x) '
-
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <subroutine walknorm calulates the distance in coordinate space for each set
-      subroutine getwalknorm(x,walknorm,npt)
-      use dim_parameter,only: qn
-      implicit none
-!     IN: variables
-      integer npt
-      double precision x(qn,npt)
-      double precision walknorm
-!     INTERNAL: variables
-      double precision xdiff(qn)
-      integer i
-
-      walknorm =0.d0
-      do i=2,npt
-         xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1)
-         walknorm = walknorm + dsqrt(sum(xdiff(1:qn)**2))
-      enddo
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <Subroutine for generating output filenames and openeing the correspondign files
-      subroutine open_plotfile(id_plot,set)
-      implicit none
-!     IN: Variables
-      integer id_plot,set
-!     INTERNAL: Variables
-      character(len=30) name    !name of output file
-
-!     define name sheme for plot files
-      if (set .lt. 10 ) then 
-            write(name,203) set 
-      else
-            write(name,202) set 
-      endif
-      
- 202  format('scan',I2,'.dat')
- 203  format('scan0',I1,'.dat')
-      !write (name,202) set
-
-c     open plotfile
-      open(id_plot,file=name)
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <Subroutine for generating output filenames and openeing the correspondign files
-      subroutine open_outfile(id_out,it_makro)
-      implicit none
-      integer id_out,it_makro
-      character(len=30) outname !name of output file
-
- 543  format('mnlfit-',i1,'.out')
- 544  format('mnlfit-',i2,'.out')
- 545  format('mnlfit-',i3,'.out')
-
-      if(it_makro.lt.10) then
-         write(outname,543) it_makro
-      else if (it_makro.lt.100) then
-         write(outname,544) it_makro
-      else if (it_makro.lt.1000) then
-         write(outname,545) it_makro
-      else
-         write(6,*)
-     >        'ERROR: No rule for Outputfile naming for MAXIT >= 1000'
-         stop
-      endif
-
-      open (id_out,file=outname)
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <Subroutine for printing the Parameterkeys for use in Input File
-!     < prints the keystring given in keys.incl and the corresponding parameters when there was atleast one parameter given in the input for the spcific key
-!     < how many parameters and spreads per line are printed can be specified with the hardcoded parameters np and nsp but they must be atleast >=2
-!     <@param id_out specifies the file in which the Parameters are Printed
-!     <@param p vector containing one set of parameter values
-!     <@param p_act vector containing the active state 0 (inactive) or 1 (active) for each parameter
-!     <@param p_spread vector containing the spreads for each parameter
-!     <@param npar lenght of the parmeter vectors (p,p_act,p_spread)
-!     <@TODO extract subroutine for printing the multiline values, would make this more readable
-      subroutine print_parameterblock(id_out,p,p_act,p_spread,npar)
-      use dim_parameter,only: pst, facspread
-      use io_parameters,only: key, parkeynum,parkeylen,llen
-      implicit none
-!     IN: Variables
-      integer  id_out,npar,p_act(npar)
-      double precision  p(npar),p_spread(npar)
-
-!     INTERNAL: variables
-!     loop index
-      integer i,k,l,t,n         !< internal variables for loops and positions in parameter vectors
-
-!     number of values per line, values must be atleast 2 set this to personal preference
-      integer, parameter :: np=5,nsp=5
-
-      character(len=llen) fmt
-
-
-!     Write header for Parameter block
- 1    format('!',200('='))
-      write(id_out,1)
-      write(id_out,'(A2,5x,A11,i3)') '! ','PARAMETER: ',npar
-      write(id_out,1)
-
-!     loop over all Parameter Keys
-      do i = 1, parkeynum
-!        save start and end of parameter block for specific key
-         k = pst(1,i)
-         l = pst(1,i)+pst(2,i)-1
-!        print only used keys with atleast one parameter
-         if(pst(2,i).gt.0) then
-            write(fmt,'("(a",I3,"'' ''i3)")') parkeylen
-            write(id_out,fmt) adjustl(key(1,i)), pst(2,i)
-
-!           write the actual parameters -> subroutine print_parameterlines()?
-            if(l-k.le.(np-1)) then
-               write(fmt,'("(a",I3,"'' ''",I3,"g24.15)")') parkeylen,np
-               write(id_out,fmt) key(2,i),(p(n), n=k,l)
-
-            else
-!              start of multi line parameter print, number of values per line specified by np
-               write(fmt,'("(a",I3,"'' ''",I3,"g24.15'' &'')")')
-     $              parkeylen,np
-               write(id_out,fmt) key(2,i),(p(n), n=k,k+(np-1))
-
-               t=k+np
-!              write continuation lines till left parameters fit on last line
-               do while(t.le.l)
-                  if(l-t.le.(np-1)) then
-                     write(fmt,'("(",I3,"x'' ''",I3,"g24.15)")')
-     $                    parkeylen,np
-                     write(id_out,fmt) (p(n), n=t, l)
-
-                  else
-                     write(fmt,'("(",I3,"x'' ''",I3,"g24.15'' &'')")')
-     $                    parkeylen,np
-                     write(id_out,fmt) (p(n), n=t, t+(np-1))
-
-                  endif
-                  t=t+np
-               enddo
-
-            endif               !-> end subroutine print_parameterlines
-
-!           write parameter active state in one line
-            write(fmt,'("(a",I3,"'' ''","50i3)")') parkeylen
-            write(id_out,fmt) key(3,i),(p_act(n),n=k,l)
-
-!           write the spreads for each parameter
-            if(l-k.le.(np-1)) then
-               write(fmt,'("(a",I3,"'' ''",I3,"g24.8)")') parkeylen,nsp
-               write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,l)
-
-            else
-!              start of multiline spread values
-               write(fmt,'("(a",I3,"'' ''",I3,"g24.8'' &'')")')
-     $              parkeylen,nsp
-               write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,k
-     >              +(np-1))
-
-               t=k+nsp
-!              write continuation lines till left spreads fit on last line
-               do while(t.le.l)
-                  if(l-t.le.(np-1)) then
-                     write(fmt,'("(",I3,"x'' ''",I3,"g24.8)")')
-     $                    parkeylen,nsp
-                     write(id_out,fmt) (p_spread(n)/facspread, n=t, l)
-                  else
-                     write(fmt,'("(",I3,"x'' ''",I3,"g24.8'' &'')")')
-     $                    parkeylen,nsp
-                     write(id_out,fmt) (p_spread(n)/facspread, n=t, t
-     >                    +(np-1))
-
-                  endif
-                  t=t+np
-               enddo
-
-            endif
-!           print empty line between diffrent parameter blocks for better readability
-            write(id_out,'(" ")')
-         endif
-
-      enddo
-
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <Subroutine for printing the current Parameters and their active state
-!     < prints only the numeric values of the parameters and does not specify the corresponding key
-!     <@param npar number of parameter
-!     <@param id_out specifies the output file
-!     <@param p,p_act parameter vectors containing the values and the activity state of parameters
-      subroutine print_parameterstate(id_out,p,p_act,npar)
-      implicit none
-
-!     IN: Variables
-      integer npar,id_out
-      double precision p(npar)
-      integer p_act(npar)
-
-!     INTERNAL: Variables
-      integer i                 !< loop control
-      integer nopt              !< number of counted active parameters
-      character(len=16) opt(npar) !< string for optimisation state
-
-!     initialize number of opt parameters and the string vector opt
-      nopt=0
-      opt = '   not opt.   '
-!     loop over all parameters and check their active state count if active and set string to opt
-      do i=1,npar
-!     Nicole: change due to value 2 of p_act
-!         if(p_act(i).eq.1) then
-         if(p_act(i).ge.1) then
-            opt(i) = '     opt.     '
-            nopt=nopt+1
-         endif
-      enddo
-!     print the Parameters and their active state within separating lines
-      write(id_out,*)''
-      write(id_out,block_line)
-      write(id_out,*) 'Parameters:'
-      write(id_out,sep_line)
-      write(id_out,'(5g14.6)') (p(i),i=1,npar)
-      write(id_out,'(5a14)')   (opt(i),i=1,npar)
-      write(id_out,sep_line)
-      write(id_out,'("No. of optimized parameters: ",i6)') nopt
-      write(id_out,block_line)
-      write(id_out,*)''
-      end subroutine
-
-!----------------------------------------------------------------------------------------------------
-!     <Subroutine for printing coordinates,refdata,modeldata,diffrence between them and the weights
-!     <@param id_out identiefies the output file
-!     <@param x vector of input pattern for each datapoint
-!     <@param y vector of expected output patterns for each datapoint
-!     <@param ymod vector of output patterns generated by the model depending on paramerters
-!     <@param wt vector of weights for each datapoint
-!     <@param qn number of input patterns
-!     <@param ntot total number of output patterns for each datapoint
-!     <@param numdatpt number of totatl datapoints
-!     <@param sets number of sets the datapoints are divided into
-!     <@param ndata vector containing the number of included datapoints for each set
-!     <@param i,j,point internal variables for loop controll and datapoint counting
-      subroutine print_data(id_out,x,y,ymod,wt)
-      use dim_parameter,only: sets,ndata,qn,ntot,numdatpt,qn_read
-      implicit none
-!     IN: Variables
-      integer id_out
-      double precision  x(qn,numdatpt)
-      double precision  y(ntot,numdatpt),ymod(ntot,numdatpt)
-      double precision  wt(ntot,numdatpt)
-
-!     INTERNAL: Variables
-      integer i,j,point
-
- 18   format(A8,i6)
- 19   format (3(A15,3x), 2x, A18  , 4x, A12)
-
-!     print seperating line and header for Data output
-      write(id_out,*) 'Printing Data Sets:'
-
-      write(id_out,19) adjustl('y(x)'),adjustl('ymod(x)'),
-     >     adjustl('y(x)-ymod(x)'),adjustl('weight'),
-     >     adjustl('x(1:qn_read) ')
-      write(id_out,sep_line)
-!     loop over all datapoints for each set and count the actual datapointnumber with point
-      point=0
-      do i=1,sets
-         write(id_out,18) 'Set: ', i
-         do j=1,ndata(i)
-            write(id_out,18) 'Point: ', j
-            point=point+1
-!           print all data for one datapoint
-            call print_datapoint(id_out,x(:,point),y(:,point),
-     >           ymod(:,point),wt(:,point))
-            write(id_out,sep_line)
-
-         enddo
-      enddo
-!     write end of data statement and two seperating lines
-      write(id_out,block_line)
-      write(id_out,*) ''
-      end subroutine
-!----------------------------------------------------------------------------------------------------
-!     <Subroutine prints a single Datapoint splits Data in nstat nci(ndiab) blocks for readability
-!     <@param id_out identiefies the output file
-!     <@param x vector of input pattern for each datapoint
-!     <@param y vector of expected output patterns for each datapoint
-!     <@param ymod vector of output patterns generated by the model depending on paramerters
-!     <@param wt vector of weights for each datapoint
-!     <@param qn number of input patterns
-!     <@param ntot total number of output patterns for each datapoint
-!     <@param i,j,k internal variables for loop controll and counting
-      subroutine print_datapoint(id_out,x,y,ymod,wt)
-      use dim_parameter,only: nstat,ndiab,nci,qn,ntot,qn_read
-      use io_parameters,only: llen
-      implicit none
-      integer id_out
-      double precision  x(qn),y(ntot),ymod(ntot),wt(ntot)
-
-      integer i,j,k
-
- 18   format(A10,i3)
- 19   format(3F18.8, 2X, F18.6, 4X,*(F12.6))
-
-!     print the nstat output patterns
-      do i=1,nstat
-         write(id_out,19)y(i),ymod(i),ymod(i)-y(i), wt(i), x(1:qn)
-      enddo
-!     loop over number (nci) of  metadata with lenght (ndiab)
-      do i=1,nci
-         write(id_out,18) 'nci: ',i
-         do j=1,ndiab
-            k=nstat + (i-1)*ndiab + j
-            write(id_out,19) y(k),ymod(k),(ymod(k)-y(k)),
-     >           wt(k), x(1:qn_read)
-         enddo
-      enddo
-
-      end subroutine
-
-      
-
-      end module write_mod