From 2c9c856bb400af0a4af4da6811a16eedb71c09be Mon Sep 17 00:00:00 2001 From: Jean Paul Nshuti Date: Thu, 18 Jun 2026 09:21:49 +0200 Subject: [PATCH] remove the model --- src/model/Phase_correction.f90 | 149 --- src/model/all.f90 | 1843 ------------------------------- src/model/ctrans.f | 58 - src/model/data_transform.f90 | 174 --- src/model/dipole_model.f90 | 411 ------- src/model/fit_genetic.f90 | 914 --------------- src/model/fit_param_dip_pes.f90 | 914 --------------- src/model/invariants_nh3.f90 | 172 --- src/model/newmodes.f | 87 -- src/model/nnadia_nh3.f | 87 -- src/model/nncoords_nh3.f90 | 91 -- src/model/pes_model.f90 | 172 --- src/model/red_invariants.f90 | 72 -- src/model/sh-scal-50.incl | 9 - src/model/shift-scale-70N.incl | 11 - src/model/write_error.f90 | 235 ---- 16 files changed, 5399 deletions(-) delete mode 100644 src/model/Phase_correction.f90 delete mode 100644 src/model/all.f90 delete mode 100644 src/model/ctrans.f delete mode 100644 src/model/data_transform.f90 delete mode 100644 src/model/dipole_model.f90 delete mode 100644 src/model/fit_genetic.f90 delete mode 100644 src/model/fit_param_dip_pes.f90 delete mode 100644 src/model/invariants_nh3.f90 delete mode 100644 src/model/newmodes.f delete mode 100644 src/model/nnadia_nh3.f delete mode 100644 src/model/nncoords_nh3.f90 delete mode 100644 src/model/pes_model.f90 delete mode 100644 src/model/red_invariants.f90 delete mode 100644 src/model/sh-scal-50.incl delete mode 100644 src/model/shift-scale-70N.incl delete mode 100644 src/model/write_error.f90 diff --git a/src/model/Phase_correction.f90 b/src/model/Phase_correction.f90 deleted file mode 100644 index 5d0a3d0..0000000 --- a/src/model/Phase_correction.f90 +++ /dev/null @@ -1,149 +0,0 @@ -module Phase_corr - ! module for correcting the phase issue in the ab inition dipole data - - use iso_fortran_env, only: dp => real64, idp => int32 - use nncommons, only: ndiab - implicit none - - !integer(idp), parameter:: ndiab = 4 - contains - - - SUBROUTINE build_cluster(m,p) - ! subroutine which check the sign of two succesive point - !and correct it if neccessary - ! take dipole mx at point i and at point i+1 - ! update dipole at i+1 - implicit none - Real(dp), intent(in):: m(ndiab,ndiab) ! at point i - Real(dp), intent(inout):: p(ndiab,ndiab) ! mx at i+1 - !Real(dp), intent(inout):: q(qn) - !Real(dp):: prod2, prod3 - - ! check first the sign of mx0(2,3) if it the same as sign of y(coord) - ! q3 is ys - ! q5 is yb - - - - - ! check if the product of the transition element - ! of two successive point there is no flip of sign - ! check if state 1 need to be flipped - !if (needs_flip(M,P,1)) then - ! write(6,*) "state 1 is flipped" - ! call flip_e_state(P,1) - !endif - - ! check if e-state (2) need to be flipped - - if (needs_flip(M,P,2)) then - !write(6,*) "state 2 is flipped" - call flip_e_state(P,2) - endif - - ! check if state (3) need to be flipped - !if (needs_flip(M,P,3)) then - ! write(6,*) 'State 3 is flipped ' - ! call flip_e_state(P,3) - !endif - - ! check for the last state - !if (needs_flip(M,P,4)) then - ! write(6,*) "state 4 is flipped" - ! call flip_e_state(P,4) - !endif - - END SUBROUTINE build_cluster - - - SUBROUTINE flip_e_state(A,O) - ! the routine which flip the sign of the e-state element - !USE dim_parameter, only: ndiab - !USE accuracy_constants, only: idp, dp - implicit none - real(dp), intent(inout):: A(ndiab,ndiab) ! dipole matrix - integer(idp), intent(in):: O ! the index of the e-state to be flipped - real(dp):: P(ndiab,ndiab) - integer(idp):: i - !logical,save :: first_call =.true. - - ! the identity matrix will be used to create the transformation matrix P - - P = 0.0d0 - do i =1,ndiab - P(i,i)=1.0d0 - enddo - ! flipped the sign of the desired e state - - P(o,o) = -1.0d0 - - ! apply P on A - ! A' = P**TAP - A=matmul(transpose(P),matmul(A,P)) - - ! write the index of electronic state flipped - - - - END SUBROUTINE - - logical function needs_flip(m, p, s) - ! Decide whether electronic state s should be flipped - implicit none - real(dp), intent(in) :: m(ndiab, ndiab), p(ndiab, ndiab) - integer(idp), intent(in) :: s - real(dp), parameter:: tol = 1.0d-9,eps_dip =1.0d-6 - - real(dp) :: m1,m2,p1,p2,m3,p3 - - needs_flip = .false. - - select case (s) - case(1) - m1 = m(1,2) ; p1 = p(1,2) - m2 = m(1,3) ; p2 = p(1,3) - m3 = m(1,4) ; p3 = p(1,4) - if ((abs(m1) < eps_dip) .and. (abs(p1) < eps_dip)) then - m1 = m3 - p1 = p3 - else if ((abs(m2) < eps_dip) .and. (abs(p2) < eps_dip)) then - m2 = m3 - p2 = p3 - endif - - case (2) - m1 = m(2,3) ; p1 = p(2,3) - m2 = m(2,4) ; p2 = P(2,4) - - case (3) - m1 = m(2,3) ; p1 = p(2,3) - m2 = m(3,4) ; p2 = p(3,4) - - case(4) - m1 = m(2,4) ; p1 = p(2,4) - m2 = m(3,4) ; p2 = p(3,4) - m3 = m(1,4) ; p3 = p(1,4) - if ((abs(m1) < eps_dip) .and. (abs(p1) < eps_dip)) then - m1 = m3 - p1 = p3 - else if ((abs(m2) < eps_dip) .and. (abs(p2) < eps_dip)) then - m2 = m3 - p2 = p3 - endif - - case default - return - end select - - if (abs(m1*p1) < tol .and. abs(m2*p2) < tol ) return - !if (abs(m2*p2) < tol ) return - - !if (m1*p1 < 0.0_dp .and. m2*p2 < 0.0_dp) then - if (m1*p1 < 0.0_dp) then - needs_flip = .true. - end if - - end function needs_flip - -end module diff --git a/src/model/all.f90 b/src/model/all.f90 deleted file mode 100644 index 07c58e9..0000000 --- a/src/model/all.f90 +++ /dev/null @@ -1,1843 +0,0 @@ -Module invariants_mod - implicit none - contains - !---------------------------------------------------- - subroutine invariants(a,xs,ys,xb,yb,inv) - implicit none - !include "nnparams.incl" - double precision, intent(in) :: a, xs, ys, xb, yb - double precision, intent(out) ::inv(3) - double precision:: invar(23) - !double precision:: invar(24) - complex(8) :: q1, q2 - LOGICAL,PARAMETER:: debg =.FALSE. - integer :: i - ! express the coordinate in complex - - q1 = dcmplx(xs, ys) - q2 = dcmplx(xb, yb) - - - ! compute the invariants - !inv( 4) = a - !inv( 4) = b**2 - - ! INVARIANTS OF KIND II - !------------------------ - - invar(1) = dreal( q1 * conjg(q1) ) ! r11 - invar(2) = dreal( q1 * conjg(q2) ) ! r12 - invar(3) = dreal( q2 * conjg(q2) ) ! r22 - invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2 - - - !INVATIANTS OF KIND III - !------------------------ - - invar(5) = dreal( q1 * q1 * q1 ) ! r111 - invar(6) = dreal( q1 * q1 * q2 ) ! r112 - invar(7) = dreal( q1 * q2 * q2 ) ! r122 - invar(8) = dreal( q2 * q2 * q2 ) ! r222 - invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2 - invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2 - invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2 - invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222 - - ! INVARIANTS OF KIND IV - !------------------------- - - invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 )) - invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 )) - invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 )) - invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 )) - - ! INVARIANTS OF KIND V - !---------------------- - - invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 )) - invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 )) - invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 )) - invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 )) - invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 )) - invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 )) - invar(23) = a - - ! the only non zero invariant for bend pure cuts - - inv(1) = invar( 1) - inv(2) = invar( 5) - inv(3) = invar( 9) - ! inv(4) = invar(1) - ! inv(5) = invar(5) - ! inv(6) = invar(9) - !inv(17)=invar(18) - !inv(19) = invar() - !inv(2) = invar(5) - !inv(3) = invar(9) - !inv(1:22)=invar(1:22) - !inv(4:8) = invar(5:9) - !inv(9)=invar(12) - - - if (debg) then - write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) - write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) - write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) - write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) - - write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES" - write(14,*)"---------------------------------------" - write(14,*)"xs =",dreal(q1), "ys=",dimag(q1) - endif - ! modify the invariants to only consider few of them - ! - !invar(13:22)=0.0d0 - end subroutine invariants - - -! subroutine invariants_ch(a, xs, ys, xb, yb, invar) -! implicit none -! double precision, intent(in) :: a, xs, ys, xb, yb -! double precision, intent(out) :: invar(25) -! -! complex(8) :: q1, q2 -! double precision :: r11, r22, r12, eta -! double precision :: I1, I2, I3, I4 -! -! ! Define complex coordinates -! q1 = dcmplx(xs, ys) -! q2 = dcmplx(xb, yb) -! -! ! Quadratic invariants -! r11 = dreal(q1*conjg(q1)) -! r22 = dreal(q2*conjg(q2)) -! r12 = dreal(q1*conjg(q2)) -! eta = dimag(q1*conjg(q2)) -! -! ! Cubic imag parts (used many times) -! I1 = dimag(q1*q1*q1) -! I2 = dimag(q1*q1*q2) -! I3 = dimag(q1*q2*q2) -! I4 = dimag(q2*q2*q2) -! ! the totally symmetric invariant 25 -! invar(25) = a -! ! Store invariants -! invar(1) = r11 -! invar(2) = r12 -! invar(3) = r22 -! invar(4) = eta*eta -! -! invar(5) = dreal(q1*q1*q1) -! invar(6) = dreal(q1*q1*q2) -! invar(7) = dreal(q1*q2*q2) -! invar(8) = dreal(q2*q2*q2) -! -! invar(9) = I1*I1 -! invar(10) = I2*I2 -! invar(11) = I3*I3 -! invar(12) = I4*I4 -! -! invar(13) = eta * I1 -! invar(14) = eta * I2 -! invar(15) = eta * I3 -! invar(16) = eta * I4 -! -! ! Pure quartic invariants -! invar(17) = r11*r11 + r22*r22 -! invar(18) = r11*r22 - r12*r12 - eta*eta -! -! ! Cubic-imag bilinear invariants (minimal independent set) -! invar(19) = I1*I2 -! invar(20) = I1*I3 -! invar(21) = I1*I4 -! invar(22) = I2*I3 -! invar(23) = I2*I4 -! -! ! Odd cubic invariant -! invar(24) = eta**3 -! -! if (debg) then -! write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) -! write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) -! write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) -! write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) -! -! write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES" -! write(14,*)"---------------------------------------" -! write(14,*)"xs =",dreal(q1), "ys=",dimag(q1) -! endif -! -! end subroutine - -end module invariants_mod - module trans_coord - implicit none - contains - - subroutine cart2mode(qq) - use invariants_mod, only: invariants - use nn_params - use nncommons - implicit none - double precision, intent(inout):: qq(maxnin) - double precision :: q(12),inv(3),qm(6) - double precision, parameter :: m_N=25526.04237518618d0 - double precision, parameter :: m_H=1837.152647439619d0 - double precision, parameter :: Ang2Bohr=1.889726124565d0 - - double precision diff(12), ref_geom(12) - double precision q_mod(12) -! transformation matrix - double precision mode(12,6) - integer i -! final internal coord vector - - - include 'newmodes.f' - - ref_geom(1:12)=0.d0 - ref_geom( 4 )=1.022871078d0 - ref_geom( 7 )=-.5d0*ref_geom(4) - !ref_geom( 8 )=sqrt(1.5d0)*ref_geom(4) - ref_geom( 8)=0.5d0*sqrt(3.0d0)*ref_geom(4) - ref_geom(10 )=-.5d0*ref_geom(4) - !ref_geom(11 )=-sqrt(1.5d0)*ref_geom(4) - ref_geom(11)=-0.5d0*sqrt(3.0d0)*ref_geom(4) - - !massN=25526.04237518618 - !massH=1837.152647439619 - ! - q(1:12)=qq(1:12) - q_mod = (1.0d0/Ang2Bohr)*q - -! difference vector - diff=q_mod-ref_geom - -! mass weighting - diff(1:3)=diff(1:3)*sqrt(m_N) - diff(4:12)=diff(4:12)*sqrt(m_H) - - qm=matmul(transpose(mode),diff) - do i =1,size(qm) - if (abs(qm(i)) .lt. 1.0d-6) qm(i)=0.0d0 - enddo - qm(3) = -qm(3) - call invariants(qm(1),qm(2),qm(3),qm(4),qm(5),inv) - qq =0.0d0 - qq(1:len_in)=inv(1:len_in) - qq(len_in+1:len_in+6)=qm(1:6) - - - end subroutine cart2mode - end module trans_coord -!**** Define coordinate transformation applied to the input before fit. -!*** -!*** -!****Conventions: -!*** -!*** ctrans: subroutine transforming a single point in coordinate space - - subroutine trans_in(pat_in) - !use ctrans_mod, only: ctrans - use trans_coord, only: cart2mode - use nn_params - use nncommons - implicit none - - !include 'nndbg.incl' - - double precision pat_in(maxnin,maxpats) - !integer ntot - - integer j,i,jj - jj=1 - do j=1,sets - !write(62,*)"Scan ",j - do i =1,ndata(j) - call cart2mode(pat_in(:,jj)) - ! FIRST ELEMENT OF PAT-IN ARE USED BY NEURON NETWORK - !write(62,'(*(ES16.8))') pat_in(:,jj) - jj = jj +1 - enddo - enddo - end - Module dip_param - IMPLICIT NONE - Integer,parameter :: np=101 - Double precision :: p(101) - integer :: pst(2,400) - ! Non zero parameters are 17 - !$omp threadprivate(p,pst) - contains - - SUBROUTINE init_dip_planar_data() - implicit none - p( 1)= -0.561433980D+02 - p( 2)= 0.000000000D+00 - p( 3)= 0.000000000D+00 - p( 4)= -0.558830090D+02 - p( 5)= 0.000000000D+00 - p( 6)= 0.000000000D+00 - p( 7)= -0.557220030D+02 - p( 8)= 0.000000000D+00 - p( 9)= 0.000000000D+00 - p( 10)= 0.000000000D+00 - p( 11)= 0.334452655D-02 - p( 12)= 0.000000000D+00 - p( 13)= 0.000000000D+00 - p( 14)= 0.000000000D+00 - p( 15)= 0.000000000D+00 - p( 16)= 0.000000000D+00 - p( 17)= 0.000000000D+00 - p( 18)= 0.000000000D+00 - p( 19)= 0.000000000D+00 - p( 20)= 0.000000000D+00 - p( 21)= -0.532590808D-03 - p( 22)= 0.000000000D+00 - p( 23)= 0.000000000D+00 - p( 24)= 0.000000000D+00 - p( 25)= 0.000000000D+00 - p( 26)= 0.000000000D+00 - p( 27)= 0.000000000D+00 - p( 28)= -0.226116045D-01 - p( 29)= 0.000000000D+00 - p( 30)= -0.363599361D-01 - p( 31)= 0.000000000D+00 - p( 32)= 0.962100323D-01 - p( 33)= 0.000000000D+00 - p( 34)= 0.111336127D-03 - p( 35)= 0.000000000D+00 - p( 36)= 0.000000000D+00 - p( 37)= -0.433685770D-04 - p( 38)= 0.000000000D+00 - p( 39)= 0.000000000D+00 - p( 40)= -0.548870572D-02 - p( 41)= 0.000000000D+00 - p( 42)= 0.000000000D+00 - p( 43)= 0.000000000D+00 - p( 44)= 0.000000000D+00 - p( 45)= 0.000000000D+00 - p( 46)= 0.000000000D+00 - p( 47)= 0.000000000D+00 - p( 48)= 0.000000000D+00 - p( 49)= 0.217598103D+00 - p( 50)= -0.113909871D-01 - p( 51)= 0.000000000D+00 - p( 52)= -0.123184224D-03 - p( 53)= 0.000000000D+00 - p( 54)= 0.000000000D+00 - p( 55)= 0.000000000D+00 - p( 56)= 0.000000000D+00 - p( 57)= 0.000000000D+00 - p( 58)= 0.000000000D+00 - p( 59)= 0.000000000D+00 - p( 60)= 0.000000000D+00 - p( 61)= 0.000000000D+00 - p( 62)= 0.000000000D+00 - p( 63)= 0.000000000D+00 - p( 64)= 0.000000000D+00 - p( 65)= 0.000000000D+00 - p( 66)= 0.000000000D+00 - p( 67)= 0.000000000D+00 - p( 68)= 0.000000000D+00 - p( 69)= 0.000000000D+00 - p( 70)= 0.000000000D+00 - p( 71)= 0.000000000D+00 - p( 72)= 0.000000000D+00 - p( 73)= 0.000000000D+00 - p( 74)= 0.000000000D+00 - p( 75)= 0.000000000D+00 - p( 76)= 0.000000000D+00 - p( 77)= 0.000000000D+00 - p( 78)= 0.000000000D+00 - p( 79)= 0.000000000D+00 - p( 80)= 0.000000000D+00 - p( 81)= 0.000000000D+00 - p( 82)= 0.000000000D+00 - p( 83)= 0.000000000D+00 - p( 84)= 0.000000000D+00 - p( 85)= 0.000000000D+00 - p( 86)= 0.000000000D+00 - p( 87)= -0.567206500D-02 - p( 88)= -0.613490432D-02 - p( 89)= 0.000000000D+00 - p( 90)= 0.303022045D-03 - p( 91)= 0.000000000D+00 - p( 92)= 0.000000000D+00 - p( 93)= 0.000000000D+00 - p( 94)= 0.000000000D+00 - p( 95)= 0.000000000D+00 - p( 96)= 0.000000000D+00 - p( 97)= 0.000000000D+00 - p( 98)= 0.000000000D+00 - p( 99)= 0.000000000D+00 - p(100)= 0.000000000D+00 - p(101)= 0.000000000D+00 - pst(1, 1)= 1 - pst(2, 1)= 3 - pst(1, 2)= 4 - pst(2, 2)= 3 - pst(1, 3)= 7 - pst(2, 3)= 4 - pst(1, 4)= 11 - pst(2, 4)= 2 - pst(1, 5)= 13 - pst(2, 5)= 3 - pst(1, 6)= 16 - pst(2, 6)= 2 - pst(1, 7)= 18 - pst(2, 7)= 3 - pst(1, 8)= 21 - pst(2, 8)= 3 - pst(1, 9)= 24 - pst(2, 9)= 3 - pst(1, 10)= 27 - pst(2, 10)= 1 - pst(1, 11)= 28 - pst(2, 11)= 2 - pst(1, 12)= 30 - pst(2, 12)= 2 - pst(1, 13)= 32 - pst(2, 13)= 2 - pst(1, 14)= 34 - pst(2, 14)= 3 - pst(1, 15)= 37 - pst(2, 15)= 3 - pst(1, 16)= 40 - pst(2, 16)= 3 - pst(1, 17)= 43 - pst(2, 17)= 2 - pst(1, 18)= 45 - pst(2, 18)= 2 - pst(1, 19)= 47 - pst(2, 19)= 2 - pst(1, 20)= 49 - pst(2, 20)= 1 - pst(1, 21)= 50 - pst(2, 21)= 2 - pst(1, 22)= 52 - pst(2, 22)= 5 - pst(1, 23)= 57 - pst(2, 23)= 10 - pst(1, 24)= 67 - pst(2, 24)= 1 - pst(1, 25)= 68 - pst(2, 25)= 2 - pst(1, 26)= 70 - pst(2, 26)= 4 - pst(1, 27)= 74 - pst(2, 27)= 8 - pst(1, 28)= 82 - pst(2, 28)= 2 - pst(1, 29)= 84 - pst(2, 29)= 3 - pst(1, 30)= 87 - pst(2, 30)= 1 - pst(1, 31)= 88 - pst(2, 31)= 2 - pst(1, 32)= 90 - pst(2, 32)= 4 - pst(1, 33)= 94 - pst(2, 33)= 8 - pst(1, 34)= 0 - pst(2, 34)= 0 - pst(1, 35)= 0 - pst(2, 35)= 0 - pst(1, 36)= 0 - pst(2, 36)= 0 - pst(1, 37)= 0 - pst(2, 37)= 0 - pst(1, 38)= 0 - pst(2, 38)= 0 - pst(1, 39)= 0 - pst(2, 39)= 0 - pst(1, 40)= 0 - pst(2, 40)= 0 - pst(1, 41)= 0 - pst(2, 41)= 0 - pst(1, 42)= 0 - pst(2, 42)= 0 - pst(1, 43)= 0 - pst(2, 43)= 0 - pst(1, 44)= 0 - pst(2, 44)= 0 - pst(1, 45)= 0 - pst(2, 45)= 0 - pst(1, 46)= 0 - pst(2, 46)= 0 - pst(1, 47)= 0 - pst(2, 47)= 0 - pst(1, 48)= 0 - pst(2, 48)= 0 - pst(1, 49)= 0 - pst(2, 49)= 0 - pst(1, 50)= 0 - pst(2, 50)= 0 - pst(1, 51)= 0 - pst(2, 51)= 0 - pst(1, 52)= 0 - pst(2, 52)= 0 - pst(1, 53)= 0 - pst(2, 53)= 0 - pst(1, 54)= 0 - pst(2, 54)= 0 - pst(1, 55)= 0 - pst(2, 55)= 0 - pst(1, 56)= 0 - pst(2, 56)= 0 - pst(1, 57)= 0 - pst(2, 57)= 0 - pst(1, 58)= 0 - pst(2, 58)= 0 - pst(1, 59)= 0 - pst(2, 59)= 0 - pst(1, 60)= 0 - pst(2, 60)= 0 - pst(1, 61)= 0 - pst(2, 61)= 0 - pst(1, 62)= 0 - pst(2, 62)= 0 - pst(1, 63)= 0 - pst(2, 63)= 0 - pst(1, 64)= 0 - pst(2, 64)= 0 - pst(1, 65)= 0 - pst(2, 65)= 0 - pst(1, 66)= 0 - pst(2, 66)= 0 - pst(1, 67)= 0 - pst(2, 67)= 0 - pst(1, 68)= 0 - pst(2, 68)= 0 - pst(1, 69)= 0 - pst(2, 69)= 0 - pst(1, 70)= 0 - pst(2, 70)= 0 - pst(1, 71)= 0 - pst(2, 71)= 0 - pst(1, 72)= 0 - pst(2, 72)= 0 - pst(1, 73)= 0 - pst(2, 73)= 0 - pst(1, 74)= 0 - pst(2, 74)= 0 - pst(1, 75)= 0 - pst(2, 75)= 0 - pst(1, 76)= 0 - pst(2, 76)= 0 - pst(1, 77)= 0 - pst(2, 77)= 0 - pst(1, 78)= 0 - pst(2, 78)= 0 - pst(1, 79)= 0 - pst(2, 79)= 0 - pst(1, 80)= 0 - pst(2, 80)= 0 - pst(1, 81)= 0 - pst(2, 81)= 0 - pst(1, 82)= 0 - pst(2, 82)= 0 - pst(1, 83)= 0 - pst(2, 83)= 0 - pst(1, 84)= 0 - pst(2, 84)= 0 - pst(1, 85)= 0 - pst(2, 85)= 0 - pst(1, 86)= 0 - pst(2, 86)= 0 - pst(1, 87)= 0 - pst(2, 87)= 0 - pst(1, 88)= 0 - pst(2, 88)= 0 - pst(1, 89)= 0 - pst(2, 89)= 0 - pst(1, 90)= 0 - pst(2, 90)= 0 - pst(1, 91)= 0 - pst(2, 91)= 0 - pst(1, 92)= 0 - pst(2, 92)= 0 - pst(1, 93)= 0 - pst(2, 93)= 0 - pst(1, 94)= 0 - pst(2, 94)= 0 - pst(1, 95)= 0 - pst(2, 95)= 0 - pst(1, 96)= 0 - pst(2, 96)= 0 - pst(1, 97)= 0 - pst(2, 97)= 0 - pst(1, 98)= 0 - pst(2, 98)= 0 - pst(1, 99)= 0 - pst(2, 99)= 0 - pst(1,100)= 0 - pst(2,100)= 0 - pst(1,101)= 0 - pst(2,101)= 0 - pst(1,102)= 0 - pst(2,102)= 0 - pst(1,103)= 0 - pst(2,103)= 0 - pst(1,104)= 0 - pst(2,104)= 0 - pst(1,105)= 0 - pst(2,105)= 0 - pst(1,106)= 0 - pst(2,106)= 0 - pst(1,107)= 0 - pst(2,107)= 0 - pst(1,108)= 0 - pst(2,108)= 0 - pst(1,109)= 0 - pst(2,109)= 0 - pst(1,110)= 0 - pst(2,110)= 0 - pst(1,111)= 0 - pst(2,111)= 0 - pst(1,112)= 0 - pst(2,112)= 0 - pst(1,113)= 0 - pst(2,113)= 0 - pst(1,114)= 0 - pst(2,114)= 0 - pst(1,115)= 0 - pst(2,115)= 0 - pst(1,116)= 0 - pst(2,116)= 0 - pst(1,117)= 0 - pst(2,117)= 0 - pst(1,118)= 0 - pst(2,118)= 0 - pst(1,119)= 0 - pst(2,119)= 0 - pst(1,120)= 0 - pst(2,120)= 0 - pst(1,121)= 0 - pst(2,121)= 0 - pst(1,122)= 0 - pst(2,122)= 0 - pst(1,123)= 0 - pst(2,123)= 0 - pst(1,124)= 0 - pst(2,124)= 0 - pst(1,125)= 0 - pst(2,125)= 0 - pst(1,126)= 0 - pst(2,126)= 0 - pst(1,127)= 0 - pst(2,127)= 0 - pst(1,128)= 0 - pst(2,128)= 0 - pst(1,129)= 0 - pst(2,129)= 0 - pst(1,130)= 0 - pst(2,130)= 0 - pst(1,131)= 0 - pst(2,131)= 0 - pst(1,132)= 0 - pst(2,132)= 0 - pst(1,133)= 0 - pst(2,133)= 0 - pst(1,134)= 0 - pst(2,134)= 0 - pst(1,135)= 0 - pst(2,135)= 0 - pst(1,136)= 0 - pst(2,136)= 0 - pst(1,137)= 0 - pst(2,137)= 0 - pst(1,138)= 0 - pst(2,138)= 0 - pst(1,139)= 0 - pst(2,139)= 0 - pst(1,140)= 0 - pst(2,140)= 0 - pst(1,141)= 0 - pst(2,141)= 0 - pst(1,142)= 0 - pst(2,142)= 0 - pst(1,143)= 0 - pst(2,143)= 0 - pst(1,144)= 0 - pst(2,144)= 0 - pst(1,145)= 0 - pst(2,145)= 0 - pst(1,146)= 0 - pst(2,146)= 0 - pst(1,147)= 0 - pst(2,147)= 0 - pst(1,148)= 0 - pst(2,148)= 0 - pst(1,149)= 0 - pst(2,149)= 0 - pst(1,150)= 0 - pst(2,150)= 0 - pst(1,151)= 0 - pst(2,151)= 0 - pst(1,152)= 0 - pst(2,152)= 0 - pst(1,153)= 0 - pst(2,153)= 0 - pst(1,154)= 0 - pst(2,154)= 0 - pst(1,155)= 0 - pst(2,155)= 0 - pst(1,156)= 0 - pst(2,156)= 0 - pst(1,157)= 0 - pst(2,157)= 0 - pst(1,158)= 0 - pst(2,158)= 0 - pst(1,159)= 0 - pst(2,159)= 0 - pst(1,160)= 0 - pst(2,160)= 0 - pst(1,161)= 0 - pst(2,161)= 0 - pst(1,162)= 0 - pst(2,162)= 0 - pst(1,163)= 0 - pst(2,163)= 0 - pst(1,164)= 0 - pst(2,164)= 0 - pst(1,165)= 0 - pst(2,165)= 0 - pst(1,166)= 0 - pst(2,166)= 0 - pst(1,167)= 0 - pst(2,167)= 0 - pst(1,168)= 0 - pst(2,168)= 0 - pst(1,169)= 0 - pst(2,169)= 0 - pst(1,170)= 0 - pst(2,170)= 0 - pst(1,171)= 0 - pst(2,171)= 0 - pst(1,172)= 0 - pst(2,172)= 0 - pst(1,173)= 0 - pst(2,173)= 0 - pst(1,174)= 0 - pst(2,174)= 0 - pst(1,175)= 0 - pst(2,175)= 0 - pst(1,176)= 0 - pst(2,176)= 0 - pst(1,177)= 0 - pst(2,177)= 0 - pst(1,178)= 0 - pst(2,178)= 0 - pst(1,179)= 0 - pst(2,179)= 0 - pst(1,180)= 0 - pst(2,180)= 0 - pst(1,181)= 0 - pst(2,181)= 0 - pst(1,182)= 0 - pst(2,182)= 0 - pst(1,183)= 0 - pst(2,183)= 0 - pst(1,184)= 0 - pst(2,184)= 0 - pst(1,185)= 0 - pst(2,185)= 0 - pst(1,186)= 0 - pst(2,186)= 0 - pst(1,187)= 0 - pst(2,187)= 0 - pst(1,188)= 0 - pst(2,188)= 0 - pst(1,189)= 0 - pst(2,189)= 0 - pst(1,190)= 0 - pst(2,190)= 0 - pst(1,191)= 0 - pst(2,191)= 0 - pst(1,192)= 0 - pst(2,192)= 0 - pst(1,193)= 0 - pst(2,193)= 0 - pst(1,194)= 0 - pst(2,194)= 0 - pst(1,195)= 0 - pst(2,195)= 0 - pst(1,196)= 0 - pst(2,196)= 0 - pst(1,197)= 0 - pst(2,197)= 0 - pst(1,198)= 0 - pst(2,198)= 0 - pst(1,199)= 0 - pst(2,199)= 0 - pst(1,200)= 0 - pst(2,200)= 0 - pst(1,201)= 0 - pst(2,201)= 0 - pst(1,202)= 0 - pst(2,202)= 0 - pst(1,203)= 0 - pst(2,203)= 0 - pst(1,204)= 0 - pst(2,204)= 0 - pst(1,205)= 0 - pst(2,205)= 0 - pst(1,206)= 0 - pst(2,206)= 0 - pst(1,207)= 0 - pst(2,207)= 0 - pst(1,208)= 0 - pst(2,208)= 0 - pst(1,209)= 0 - pst(2,209)= 0 - pst(1,210)= 0 - pst(2,210)= 0 - pst(1,211)= 0 - pst(2,211)= 0 - pst(1,212)= 0 - pst(2,212)= 0 - pst(1,213)= 0 - pst(2,213)= 0 - pst(1,214)= 0 - pst(2,214)= 0 - pst(1,215)= 0 - pst(2,215)= 0 - pst(1,216)= 0 - pst(2,216)= 0 - pst(1,217)= 0 - pst(2,217)= 0 - pst(1,218)= 0 - pst(2,218)= 0 - pst(1,219)= 0 - pst(2,219)= 0 - pst(1,220)= 0 - pst(2,220)= 0 - pst(1,221)= 0 - pst(2,221)= 0 - pst(1,222)= 0 - pst(2,222)= 0 - pst(1,223)= 0 - pst(2,223)= 0 - pst(1,224)= 0 - pst(2,224)= 0 - pst(1,225)= 0 - pst(2,225)= 0 - pst(1,226)= 0 - pst(2,226)= 0 - pst(1,227)= 0 - pst(2,227)= 0 - pst(1,228)= 0 - pst(2,228)= 0 - pst(1,229)= 0 - pst(2,229)= 0 - pst(1,230)= 0 - pst(2,230)= 0 - pst(1,231)= 0 - pst(2,231)= 0 - pst(1,232)= 0 - pst(2,232)= 0 - pst(1,233)= 0 - pst(2,233)= 0 - pst(1,234)= 0 - pst(2,234)= 0 - pst(1,235)= 0 - pst(2,235)= 0 - pst(1,236)= 0 - pst(2,236)= 0 - pst(1,237)= 0 - pst(2,237)= 0 - pst(1,238)= 0 - pst(2,238)= 0 - pst(1,239)= 0 - pst(2,239)= 0 - pst(1,240)= 0 - pst(2,240)= 0 - pst(1,241)= 0 - pst(2,241)= 0 - pst(1,242)= 0 - pst(2,242)= 0 - pst(1,243)= 0 - pst(2,243)= 0 - pst(1,244)= 0 - pst(2,244)= 0 - pst(1,245)= 0 - pst(2,245)= 0 - pst(1,246)= 0 - pst(2,246)= 0 - pst(1,247)= 0 - pst(2,247)= 0 - pst(1,248)= 0 - pst(2,248)= 0 - pst(1,249)= 0 - pst(2,249)= 0 - pst(1,250)= 0 - pst(2,250)= 0 - pst(1,251)= 0 - pst(2,251)= 0 - pst(1,252)= 0 - pst(2,252)= 0 - pst(1,253)= 0 - pst(2,253)= 0 - pst(1,254)= 0 - pst(2,254)= 0 - pst(1,255)= 0 - pst(2,255)= 0 - pst(1,256)= 0 - pst(2,256)= 0 - pst(1,257)= 0 - pst(2,257)= 0 - pst(1,258)= 0 - pst(2,258)= 0 - pst(1,259)= 0 - pst(2,259)= 0 - pst(1,260)= 0 - pst(2,260)= 0 - pst(1,261)= 0 - pst(2,261)= 0 - pst(1,262)= 0 - pst(2,262)= 0 - pst(1,263)= 0 - pst(2,263)= 0 - pst(1,264)= 0 - pst(2,264)= 0 - pst(1,265)= 0 - pst(2,265)= 0 - pst(1,266)= 0 - pst(2,266)= 0 - pst(1,267)= 0 - pst(2,267)= 0 - pst(1,268)= 0 - pst(2,268)= 0 - pst(1,269)= 0 - pst(2,269)= 0 - pst(1,270)= 0 - pst(2,270)= 0 - pst(1,271)= 0 - pst(2,271)= 0 - pst(1,272)= 0 - pst(2,272)= 0 - pst(1,273)= 0 - pst(2,273)= 0 - pst(1,274)= 0 - pst(2,274)= 0 - pst(1,275)= 0 - pst(2,275)= 0 - pst(1,276)= 0 - pst(2,276)= 0 - pst(1,277)= 0 - pst(2,277)= 0 - pst(1,278)= 0 - pst(2,278)= 0 - pst(1,279)= 0 - pst(2,279)= 0 - pst(1,280)= 0 - pst(2,280)= 0 - pst(1,281)= 0 - pst(2,281)= 0 - pst(1,282)= 0 - pst(2,282)= 0 - pst(1,283)= 0 - pst(2,283)= 0 - pst(1,284)= 0 - pst(2,284)= 0 - pst(1,285)= 0 - pst(2,285)= 0 - pst(1,286)= 0 - pst(2,286)= 0 - pst(1,287)= 0 - pst(2,287)= 0 - pst(1,288)= 0 - pst(2,288)= 0 - pst(1,289)= 0 - pst(2,289)= 0 - pst(1,290)= 0 - pst(2,290)= 0 - pst(1,291)= 0 - pst(2,291)= 0 - pst(1,292)= 0 - pst(2,292)= 0 - pst(1,293)= 0 - pst(2,293)= 0 - pst(1,294)= 0 - pst(2,294)= 0 - pst(1,295)= 0 - pst(2,295)= 0 - pst(1,296)= 0 - pst(2,296)= 0 - pst(1,297)= 0 - pst(2,297)= 0 - pst(1,298)= 0 - pst(2,298)= 0 - pst(1,299)= 0 - pst(2,299)= 0 - pst(1,300)= 0 - pst(2,300)= 0 - pst(1,301)= 0 - pst(2,301)= 0 - pst(1,302)= 0 - pst(2,302)= 0 - pst(1,303)= 0 - pst(2,303)= 0 - pst(1,304)= 0 - pst(2,304)= 0 - pst(1,305)= 0 - pst(2,305)= 0 - pst(1,306)= 0 - pst(2,306)= 0 - pst(1,307)= 0 - pst(2,307)= 0 - pst(1,308)= 0 - pst(2,308)= 0 - pst(1,309)= 0 - pst(2,309)= 0 - pst(1,310)= 0 - pst(2,310)= 0 - pst(1,311)= 0 - pst(2,311)= 0 - pst(1,312)= 0 - pst(2,312)= 0 - pst(1,313)= 0 - pst(2,313)= 0 - pst(1,314)= 0 - pst(2,314)= 0 - pst(1,315)= 0 - pst(2,315)= 0 - pst(1,316)= 0 - pst(2,316)= 0 - pst(1,317)= 0 - pst(2,317)= 0 - pst(1,318)= 0 - pst(2,318)= 0 - pst(1,319)= 0 - pst(2,319)= 0 - pst(1,320)= 0 - pst(2,320)= 0 - pst(1,321)= 0 - pst(2,321)= 0 - pst(1,322)= 0 - pst(2,322)= 0 - pst(1,323)= 0 - pst(2,323)= 0 - pst(1,324)= 0 - pst(2,324)= 0 - pst(1,325)= 0 - pst(2,325)= 0 - pst(1,326)= 0 - pst(2,326)= 0 - pst(1,327)= 0 - pst(2,327)= 0 - pst(1,328)= 0 - pst(2,328)= 0 - pst(1,329)= 0 - pst(2,329)= 0 - pst(1,330)= 0 - pst(2,330)= 0 - pst(1,331)= 0 - pst(2,331)= 0 - pst(1,332)= 0 - pst(2,332)= 0 - pst(1,333)= 0 - pst(2,333)= 0 - pst(1,334)= 0 - pst(2,334)= 0 - pst(1,335)= 0 - pst(2,335)= 0 - pst(1,336)= 0 - pst(2,336)= 0 - pst(1,337)= 0 - pst(2,337)= 0 - pst(1,338)= 0 - pst(2,338)= 0 - pst(1,339)= 0 - pst(2,339)= 0 - pst(1,340)= 0 - pst(2,340)= 0 - pst(1,341)= 0 - pst(2,341)= 0 - pst(1,342)= 0 - pst(2,342)= 0 - pst(1,343)= 0 - pst(2,343)= 0 - pst(1,344)= 0 - pst(2,344)= 0 - pst(1,345)= 0 - pst(2,345)= 0 - pst(1,346)= 0 - pst(2,346)= 0 - pst(1,347)= 0 - pst(2,347)= 0 - pst(1,348)= 0 - pst(2,348)= 0 - pst(1,349)= 0 - pst(2,349)= 0 - pst(1,350)= 0 - pst(2,350)= 0 - pst(1,351)= 0 - pst(2,351)= 0 - pst(1,352)= 0 - pst(2,352)= 0 - pst(1,353)= 0 - pst(2,353)= 0 - pst(1,354)= 0 - pst(2,354)= 0 - pst(1,355)= 0 - pst(2,355)= 0 - pst(1,356)= 0 - pst(2,356)= 0 - pst(1,357)= 0 - pst(2,357)= 0 - pst(1,358)= 0 - pst(2,358)= 0 - pst(1,359)= 0 - pst(2,359)= 0 - pst(1,360)= 0 - pst(2,360)= 0 - pst(1,361)= 0 - pst(2,361)= 0 - pst(1,362)= 0 - pst(2,362)= 0 - pst(1,363)= 0 - pst(2,363)= 0 - pst(1,364)= 0 - pst(2,364)= 0 - pst(1,365)= 0 - pst(2,365)= 0 - pst(1,366)= 0 - pst(2,366)= 0 - pst(1,367)= 0 - pst(2,367)= 0 - pst(1,368)= 0 - pst(2,368)= 0 - pst(1,369)= 0 - pst(2,369)= 0 - pst(1,370)= 0 - pst(2,370)= 0 - pst(1,371)= 0 - pst(2,371)= 0 - pst(1,372)= 0 - pst(2,372)= 0 - pst(1,373)= 0 - pst(2,373)= 0 - pst(1,374)= 0 - pst(2,374)= 0 - pst(1,375)= 0 - pst(2,375)= 0 - pst(1,376)= 0 - pst(2,376)= 0 - pst(1,377)= 0 - pst(2,377)= 0 - pst(1,378)= 0 - pst(2,378)= 0 - pst(1,379)= 0 - pst(2,379)= 0 - pst(1,380)= 0 - pst(2,380)= 0 - pst(1,381)= 0 - pst(2,381)= 0 - pst(1,382)= 0 - pst(2,382)= 0 - pst(1,383)= 0 - pst(2,383)= 0 - pst(1,384)= 0 - pst(2,384)= 0 - pst(1,385)= 0 - pst(2,385)= 0 - pst(1,386)= 0 - pst(2,386)= 0 - pst(1,387)= 0 - pst(2,387)= 0 - pst(1,388)= 0 - pst(2,388)= 0 - pst(1,389)= 0 - pst(2,389)= 0 - pst(1,390)= 0 - pst(2,390)= 0 - pst(1,391)= 0 - pst(2,391)= 0 - pst(1,392)= 0 - pst(2,392)= 0 - pst(1,393)= 0 - pst(2,393)= 0 - pst(1,394)= 0 - pst(2,394)= 0 - pst(1,395)= 0 - pst(2,395)= 0 - pst(1,396)= 0 - pst(2,396)= 0 - pst(1,397)= 0 - pst(2,397)= 0 - pst(1,398)= 0 - pst(2,398)= 0 - pst(1,399)= 0 - pst(2,399)= 0 - pst(1,400)= 0 - pst(2,400)= 0 -End SUBROUTINE init_dip_planar_data -End Module dip_param -module pesmod - !use dim_parameter,only:qn,ndiab,pst - use iso_fortran_env, only: idp => int32, dp => real64 - use dip_param - - use nn_params - use nncommons - implicit none - - integer(idp),parameter:: ndiab=4 - private - public :: potential_nh3_low - contains - - - subroutine potential_nh3_low(e,q,nn_out) - implicit none - real(dp),intent(in)::q(maxnin) - real(dp),intent(out)::e(:,:) - real(dp),intent(inout):: nn_out(maxnout) - !logical,intent(in):: sym - !integer(idp),intent(in)::key - integer(idp) id,i, ii - real(dp) tmp,xs,xb,ys,yb,a,b,ss,sb - logical, parameter:: switch = .false. ! parameter to switch when you want to know Ex and Ey - logical, parameter:: sym = .false. ! to include r**2 in the Genetic fitting - real(dp),parameter:: tol = 1.0d-9 - real(dp),dimension(17):: scal, shift - - xs=q(len_in+2) - ys=q(len_in+3) - xb=q(len_in+4) - yb=q(len_in+5) - a=q(len_in+1) - b=q(len_in+6) - - ss= xs**2+ys**2 - sb= xb**2+yb**2 - call init_dip_planar_data() - scal = 1.0d-3 - shift = 1.0d0 - ii = 1 - do i =1 ,27 ! pes parameter key index - if( abs(p(i)) .gt. tol ) then - p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii)) - ii = ii + 1 - - endif - enddo - !write(6,*)" PES non zero", ii-1 - - - e = 0.0_dp - id = 1 ! 1 - ! v -term - ! oth order only - ! A2" - e(1,1) = e(1,1) + p(pst(1,id)) !- p(pst(1,1)) - if (sym) then - e(1,1)=e(1,1) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*sb - endif - id = id +1 !2 - ! E block - e(2,2) = e(2,2)+ p(pst(1,id))! - p(pst(1,1)) - e(3,3) = e(3,3)+ p(pst(1,id))! - p(pst(1,1)) - if (sym) then - e(2,2)=e(2,2)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 - e(3,3)=e(3,3)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 - endif - - id = id +1 ! 3 - - !A1' - e(4,4) = e(4,4) + p(pst(1,id))! - p(pst(1,1)) - if (sym) then - e(4,4)=e(4,4) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 & - + p(pst(1,id)+3)*ss**3 - endif - - ! W and z of JT block - - ! first order - id = id +1 !4 - - 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 ! 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) - e(2,3) = e(2,3) - p(pst(1,id))*(2*xs*ys) & - - p(pst(1,id)+1)*(2*xb*yb) & - - p(pst(1,id)+2)*(xs*yb+xb*ys) - - ! the coupling of A2' and E' - ! order 1 - id = id +1 !6 - - e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*yb) - e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) - ! order 2 - id =id +1 ! 7 - e(1,2) = e(1,2) +b*(p(pst(1,id))*(xs**2-ys**2) & - + p(pst(1,id)+1)*(xb**2-yb**2) & - + p(pst(1,id)+2)*(xs*xb -ys*yb)) - e(1,3) = e(1,3) +b*(p(pst(1,id))*(2*xs*ys) & - + p(pst(1,id)+1)*(2*xb*yb) & - + p(pst(1,id)+2)*(xs*yb + xb*ys)) - - ! the coupling of A1' and E' - id = id +1 ! 8 - e(2,4) = e(2,4) + p(pst(1,id))*xs + p(pst(1,id)+1)*xb - e(3,4) = e(3,4) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb - if (sym) then - e(2,4) = e(2,4)+ p(pst(1,id)+2)*xs*ss - e(3,4) = e(3,4)- p(pst(1,id)+2)*ys*ss - endif - - ! order 2 - id = id +1 ! 9 - e(2,4) = e(2,4) + p(pst(1,id))*(xs**2-ys**2) & - + p(pst(1,id)+1)*(xb**2-yb**2) & - + p(pst(1,id)+2)*(xs*xb -ys*yb) - e(3,4) = e(3,4) + p(pst(1,id))*(2*xs*ys) & - + p(pst(1,id)+1)*(2*xb*yb) & - + p(pst(1,id)+2)*(xs*yb + xb*ys) - ! the couplign between A - id = id +1 ! 10 - e(1,4) = e(1,4) + b*p(pst(1,id)) - !write(*,*) 'Id = ',id - - call copy_2_lower_triangle(e) - if (switch) then - write( 6,*)" EX AND EY ARE SWITCHED" - tmp = e(2,2) - e(2,2) = e(3,3) - e(3,3) = tmp - - ! the coupling - - tmp = e(2,4) - e(2,4) = e(3,4) - e(3,4) = tmp - e(4,2) = e(2,4) - e(4,3) = e(3,4) - endif - end subroutine potential_nh3_low - - - subroutine copy_2_lower_triangle(mat) - real(dp), intent(inout) :: mat(:, :) - integer :: m, n - ! write lower triangle of matrix symmetrical - do n = 2, size(mat, 1) - do m = 1, n - 1 - mat(n, m) = mat(m, n) - end do - end do - end subroutine copy_2_lower_triangle - - -end module -module diabmodel - use dip_param - use iso_fortran_env, only:dp => real64,idp => int32 - use nn_params - use nncommons - implicit none - - !integer(idp),parameter:: ndiab=4 - !logical :: debug=.false. - private - public:: diab_x - contains - - - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine diab_x(e,q,nn_out) - real(dp),intent(in) :: q(maxnin) - real(dp),intent(inout) :: nn_out(maxnout) - real(dp),intent(out) :: e(:,:) - integer(idp) id,i,j, ii - real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) - real(dp), dimension(17) :: shift,scal - real(dp), parameter:: tol = 1.0d-9 - - xs=q(len_in+2) - ys=q(len_in+3) - xb=q(len_in+4) - yb=q(len_in+5) - a=q(len_in+1) - b=q(len_in+6) - - id = 11 - - call init_dip_planar_data() - scal = 1.0d0 - !scal(19:20) = 1.0d-2 - - shift =1.0d0 - ii = 5 - do i = 28 ,np ! pes parameter key index - if( abs(p(i)) .gt. tol ) then - p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii)) - ii = ii + 1 - endif - enddo - !write(6,*)" Dipole non-zero", ii-1 - - ss=xs**2+ys**2 ! totaly symmetric term - sb=xb**2+yb**2 - - v3_vec( 1) = xs*(xs**2-3*ys**2) - v3_vec( 2) = xb*(xb**2-3*yb**2) - v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys - v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb - v3_vec( 5) = ys*(3*xs**2-ys**2) - v3_vec( 6) = yb*(3*xb**2-yb**2) - v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys - v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb - - - - e=0.0d0 - - - - !id=key !1 - ! V-term - ! order 1 - e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - id=id+1 !2 - e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - id=id+1 !3 - e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - ! order 2 - id=id+1 !4 - e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb) - id=id+1 !5 - e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +& - p(pst(1,id)+2)*(xs*xb-ys*yb) - e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +& - p(pst(1,id)+2)*(xs*xb-ys*yb) - id=id+1 !6 - e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) + & - p(pst(1,id)+2)*(xs*xb-ys*yb) - ! order 3 - id=id+1 !7 - e(1,1)=e(1,1)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - id=id+1 !8 - e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - id=id+1 !9 - e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - - ! JAHN TELLER COUPLING W AND Z - ! order 0 - id=id+1 !10 - e(2,2)=e(2,2)+p(pst(1,id)) - e(3,3)=e(3,3)-p(pst(1,id)) - ! order 1 - id=id+1 !11 - e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - e(3,3)=e(3,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb - e(2,3)=e(2,3)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb - ! order 2 - id=id+1 !12 - e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb - e(3,3)=e(3,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb) - e(2,3)=e(2,3)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+ & - p(pst(1,id)+2)*(xs*yb+xb*ys) - ! order 3 - - id=id+1 !13 - do i=1,4 - j=i-1 - e(2,2)=e(2,2)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) - e(3,3)=e(3,3)-(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) - e(2,3)=e(2,3)+(-p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i+4) - enddo - - e(2,2)=e(2,2)+p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb - e(3,3)=e(3,3)-(p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb) - e(2,3)=e(2,3)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb - ! PSEUDO JAHN TELLER - - ! A2 ground state coupled with E - ! ################################################### - ! ################################################### - - ! order 0 - id=id+1 !14 - e(1,2)=e(1,2)+b*p(pst(1,id)) - - ! order 1 - id=id+1 !15 - e(1,2)=e(1,2)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) - e(1,3)=e(1,3)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) - ! order 2 - id=id+1 !16 - e(1,2)=e(1,2)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb) + p(pst(1,id)+3)*(xs**2+ys**2)) - e(1,3)=e(1,3)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& - +p(pst(1,id)+2)*(xs*yb+xb*ys)) - ! order 3 - id =id+1 ! 17 - - do i=1,4 - e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) - e(1,3)=e(1,3)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) - enddo - - - - !! THE COUPLING OF A2 WITH A1 - !#################################################### - !#################################################### - ! order 1 - id=id+1 !18 - e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) - id=id+1 !19 - e(1,4)=e(1,4)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb)) - - - !!! THE COUPLING OF A1 WITH E - !!#################################################### - !#################################################### - ! order 0 - - id=id+1 !20 - e(2,4)=e(2,4)+p(pst(1,id)) - - ! order 1 - id=id+1 !21 - e(2,4)=e(2,4)+p(pst(1,id))*xs + p(pst(1,id)+1)*xb - e(3,4)=e(3,4)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb - - ! order 2 - id=id+1 !22 - e(2,4)=e(2,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb) +p(pst(1,id)+3)*(xs**2+ys**2) - e(3,4)=e(3,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & - -p(pst(1,id)+2)*(xs*yb+xb*ys) - ! order 3 - id=id+1 !23 - do i=1,4 - e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) - e(3,4)=e(3,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) - enddo - - !! End of the model - - e(2,1)=e(1,2) - e(3,1)=e(1,3) - e(3,2)=e(2,3) - e(4,1)=e(1,4) - e(4,2)=e(2,4) - e(4,3)=e(3,4) - end subroutine diab_x - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! THE Y COMPONENT OF DIPOLE - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!subroutine diab_y(e,q,nn_out) -! real(dp),intent(in) :: q(maxnin) -! real(dp),intent(inout) :: nn_out(maxnout) -! real(dp),intent(out) :: e(:,:) -! integer(idp) id,i,j, ii -! real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) -! real(dp), dimension(23) :: shift,scal -! real(dp), parameter:: tol = 1.0d-9 -! -! xs=q(5) -! ys=q(6) -! xb=q(7) -! yb=q(8) -! a=q(4) -! b=q(9) -! -! id = 11 -! -! call init_dip_planar_data() -! scal = 1.0d-3 -! shift =1.0d0 -! ii = 1 -! do i = 11 ,np ! pes parameter key index -! if( abs(p(i)) .gt. tol ) then -! p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii)) -! ii = ii + 1 -! endif -! enddo -! -! id = 11 -! -! ss=xs**2+ys**2 ! totaly symmetric term -! sb=xb**2+yb**2 -! -! v3_vec( 1) = xs*(xs**2-3*ys**2) -! v3_vec( 2) = xb*(xb**2-3*yb**2) -! v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys -! v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb -! v3_vec( 5) = ys*(3*xs**2-ys**2) -! v3_vec( 6) = yb*(3*xb**2-yb**2) -! v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys -! v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb -! -! -! -! -! -! e=0.0d0 -! ! V-term -! !id=key !1 -! e(1,1)=e(1,1)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! -! id=id+1 !2 -! e(2,2)=e(2,2)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! id=id+1 !3 -! e(4,4)=e(4,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! ! order 2 -! id=id+1 !4 -! e(1,1)=e(1,1)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! id=id+1 !5 -! e(2,2)=e(2,2)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! -! -! e(3,3)=e(3,3)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! id=id+1 !6 -! e(4,4)=e(4,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! ! order 3 -! id=id+1 !7 -! e(1,1)=e(1,1)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! id=id+1 !8 -! e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! id=id+1 !9 -! e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! -! ! V- term + totally symmetric coord a -! -! ! JAHN TELLER COUPLING TERM -! ! order 0 -! id=id+1 !10 -! e(2,3)=e(2,3)+p(pst(1,id)) -! ! order 1 -! -! id=id+1 !11 -! e(2,2)=e(2,2)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb -! e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! e(2,3)=e(2,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb -! !id=id+1 !12 -! ! order 2 -! id=id+1 !12 -! e(2,2)=e(2,2)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+p(pst(1,id)+2)*(xs*yb+xb*ys) -! e(3,3)=e(3,3)-p(pst(1,id))*2*xs*ys-p(pst(1,id)+1)*2*xb*yb-p(pst(1,id)+2)*(xs*yb+xb*ys) -! e(2,3)=e(2,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)) & -! -p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb -! ! order 3 -! id=id+1 !13 -! do i=1,4 -! j=i-1 -! e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) -! e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) -! e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) -! enddo -! e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb -! e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb -! e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb -! ! PSEUDO JAHN TELLER -! ! ORDER 0 -! ! THE COUPLING OF A2 GROUND STATE WITH E -! ! ################################################### -! ! ################################################### -! ! order 0 -! id=id+1 !14 -! e(1,3)=e(1,3)-b*(p(pst(1,id))) -! ! order 1 -! id=id+1 !15 -! e(1,2)=e(1,2)-b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) -! e(1,3)=e(1,3)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) -! -! ! order 2 -! id=id+1 !16 -! e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& -! +p(pst(1,id)+2)*(xs*yb+xb*ys)) -! e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& -! +p(pst(1,id)+2)*(xs*xb-ys*yb) - p(pst(1,id)+3)*(xs**2+ys**2)) -! -! ! order 3 -! id = id+1 ! 17 -! do i=1,4 -! e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) -! e(1,3)=e(1,3)-b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) -! enddo -! -! -! ! THE COUPLING OF A2 WITH A1 -! !#################################################### -! !#################################################### -! ! order 1 -! id=id+1 !17 -! e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) -! ! order 2 -! id=id+1 !18 -! e(1,4)=e(1,4)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& -! +p(pst(1,id)+2)*(xs*yb+xb*ys)) -! -! -! ! THE COUPLING OF A1 WITH E -! !#################################################### -! !#################################################### -! ! order 0 -! id=id+1 !19 -! e(3,4)=e(3,4)-p(pst(1,id)) -! ! order 1 -! id=id+1 !20 -! e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb -! e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb -! ! order 2 -! id=id+1 !21 -! e(2,4)=e(2,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb) & -! +p(pst(1,id)+2)*(xs*yb+xb*ys) -! e(3,4)=e(3,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & -! +p(pst(1,id)+2)*(xs*xb-ys*yb) - p(pst(1,id)+3)*(xs**2+ys**2) -! -! id =id+1 ! 23 -! ! order 3 -! do i=1,4 -! e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) -! e(3,4)=e(3,4)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) -! enddo -! -! ! end of the model -! e(2,1)=e(1,2) -! e(3,1)=e(1,3) -! e(3,2)=e(2,3) -! e(4,1)=e(1,4) -! e(4,2)=e(2,4) -! e(4,3)=e(3,4) -! ! change Ex and Ey -!end subroutine diab_y - -! potential nh3 -! parth - - -!subroutine copy_2_lower_triangle(mat) -! real(dp), intent(inout) :: mat(:, :) -! integer :: m, n -! write lower triangle of matrix symmetrical -! do n = 2, size(mat, 1) -! do m = 1, n - 1 -! mat(n, m) = mat(m, n) -! end do -! end do -! end subroutine copy_2_lower_triangle - -end module diabmodel - subroutine nninit_mod() - use nn_params - implicit none - - end subroutine - - - subroutine nnadia(coords,coeffs,adiaoutp) - USE pesmod, only: potential_nh3_low - Use diabmodel, only: diab_x - - use nn_params - use nncommons - implicit none - ! returns THE DIABATIC MATRIX UPPER MATRIX OF BOTH X AND Y COMPONENT OF DIPOLE -! -! coords: vector containing symmetry-adapted coordinates. -! coeffs: vector containing coeffs of diab dipole matrix - -! adiaoutp: upper diab matrix of dipole -! x component come first and then y component -! -! dmat_x & dmat_y: analytic model of dipole, x and y -! -! matdim: dimension of hamiltoniam matrix -! mkeigen: integer: no eigenvectors if =0 -! eigenvs: dummy storage for eigenvectors - - !include 'JTmod.incl' - integer ndiab - parameter (ndiab=4) - DOUBLE PRECISION,intent(in):: coords(maxnin) - DOUBLE PRECISION,INTENT(INOUT)::coeffs(maxnout) - DOUBLE PRECISION,INTENT(OUT):: adiaoutp(maxpout) - DOUBLE PRECISION, DIMENSION(ndiab,ndiab):: v,U,ex,temp - double precision, dimension(ndiab):: e - Double precision, allocatable:: mat(:,:) - integer i, ij, j, max_row - - ! variable for lapack - integer,parameter :: lwork = 1000 - double precision work(lwork) - integer info - - call potential_nh3_low(V,coords,coeffs) - - ! diagonalize the pes - u=0.0d0 - allocate(mat,source=V) - call DSYEV('V','U',ndiab,mat,ndiab,e,work,lwork, info) - - if (info .ne. 0) then - write(6,*)"The diagonalizing routine failed", info - stop - end if - U(:,:)= mat(:,:) - deallocate(mat) - do i =1,ndiab - adiaoutp(i)=e(i) - enddo - - ! Fix the phase of U matrix - ! fix phas of U - do i =1,ndiab - max_row=maxloc(abs(U(:,i)),1) - if(U(max_row,i) .lt. 0) then - U(:,i) = -1.0d0*U(:,i) - endif - enddo - - call diab_x(ex,coords,coeffs) ! ex is the dipole x component - - ! transforrm the dipole - - temp = matmul(transpose(U),ex) - ex = matmul(temp,U) - ij = 4 - do i =1,ndiab - do j = i,ndiab - ij = ij+1 - adiaoutp(ij) = ex(i,j) - enddo - enddo - - END SUBROUTINE - - - diff --git a/src/model/ctrans.f b/src/model/ctrans.f deleted file mode 100644 index 0971f33..0000000 --- a/src/model/ctrans.f +++ /dev/null @@ -1,58 +0,0 @@ -! Author: jnshuti -! Created: 2025-10-24 16:44:03 -! Last modified: 2026-01-12 11:20:55 jnshuti - subroutine cart2mode(qq) - use invariants_mod, only: invariants - use nn_params - use nncommons - implicit none - double precision, intent(inout):: qq(maxnin) - double precision :: q(12),inv(4) - double precision :: qm(6) - - double precision, parameter :: m_N=25526.04237518618d0 - double precision, parameter :: m_H=1837.152647439619d0 - double precision, parameter :: Ang2Bohr=1.889726124565d0 - - double precision diff(12), ref_geom(12) - double precision q_mod(12) -! transformation matrix - double precision mode(12,6) -! final internal coord vector - - - include 'newmodes.f' - - ref_geom(1:12)=0.d0 - ref_geom( 4 )=1.022871078d0 - ref_geom( 7 )=-.5d0*ref_geom(4) - !ref_geom( 8 )=sqrt(1.5d0)*ref_geom(4) - ref_geom( 8)=0.5d0*sqrt(3.0d0)*ref_geom(4) - ref_geom(10 )=-.5d0*ref_geom(4) - !ref_geom(11 )=-sqrt(1.5d0)*ref_geom(4) - ref_geom(11)=-0.5d0*sqrt(3.0d0)*ref_geom(4) - - !massN=25526.04237518618 - !massH=1837.152647439619 - ! - q(1:12)=qq(1:12) - q_mod = (1.0d0/Ang2Bohr)*q - -! difference vector - diff=q_mod-ref_geom - -! mass weighting - diff(1:3)=diff(1:3)*sqrt(m_N) - diff(4:12)=diff(4:12)*sqrt(m_H) - - qm=matmul(transpose(mode),diff) - do i =1,6 - if ( abs(qm(i)) .lt. 1.0d-6) qm(i) = 0.0d0 - enddo - qm(3) = -qm(3) - call invariants(qm(1),qm(2),qm(3),qm(4),qm(5),inv) - qq(1:len_in)=inv(1:len_in) - qq(len_in+1:len_in+5)=qm(2:6) - - - end subroutine cart2mode diff --git a/src/model/data_transform.f90 b/src/model/data_transform.f90 deleted file mode 100644 index 4678c3e..0000000 --- a/src/model/data_transform.f90 +++ /dev/null @@ -1,174 +0,0 @@ -! Author: jnshuti -! Created: 2025-10-28 10:05:45 -! Last modified: 2026-02-04 10:18:40 jnshuti -! Subroutine which transform the data in any desired format -module data_transf_mod - use iso_fortran_env, only: dp => real64, idp => int32 - use nncommons, only: sets, ndata, inp_out,ndiab - use nn_params, only: maxpout,maxpats - use Phase_corr - implicit none - !integer(idp),parameter:: ndiab=4 - private - public :: data_transform - contains - - SUBROUTINE data_transform(pat_out,npat) - ! IN: pat_out: the output patterns - ! x1 the input pattern - ! the aim it to transform pat_out in any disired form - implicit none - - real(dp),intent(inout):: pat_out(maxpout,maxpats) - !real(dp),intent(inout):: ref_out(maxpout,maxpats) - integer(idp), intent(in):: npat ! the total number of point - ! internal varibles - integer(idp) i,j,total ,ii,m,n - integer(idp), dimension(sets):: flipped_scan2, flipped_scan3 - integer(idp), dimension(sets):: flipped_scan4 - Real(dp),allocatable:: mux(:,:,:) - integer(idp):: scn_strt, scn_end,iref - - logical, parameter:: flip =.true. - - !write(6,*)"The total numbe of points used in data transform =",npat - - - ! allocatable varables - allocate(mux(ndiab,ndiab,npat)) - - ! define the scann need to be flipped - ! list of all the scan need to flip the sign depending on the result of the molpro calculation - flipped_scan2 = 0; flipped_scan3 = 0; flipped_scan4 =0 - flipped_scan3 = [1,3,4,5,6,7,8,9,10,11, & - 12,13,14,15,16,17,18,19,20, & - 21,22,24,25,54,55,56, & - 58,59,60,61,62,63] - flipped_scan2 = [27,28,29,30,31,32,33,34,35, & - 36,37,38,39,40,41,42,43,44, & - 45,46,47,48,49,50,51,52,53, & - 65,67,68,69,70,71, & - 72,73,74] - flipped_scan4 = [23,64] - !flipped_scan2 = [77,80,118,120] - !flipped_scan3 = [56,111] - - - - - - - - - total=0 - - do i =1,sets - do j=1,ndata(i) - total = total+1 - !pat_out(1:inp_out,total)=pat_out(1:inp_out,total)!-ref_en - ! transform the y to matrix form - - call Y2mat(pat_out(5:inp_out,total),Mux(:,:,total)) - if (flip) then - if (j .gt. 2) call build_cluster(Mux(:,:,total-1),Mux(:,:,total)) - - endif - - enddo - enddo - - - - ! flipe the required scan - if (flip) then - scn_end = 0 - iref=0 - do i =1,sets - if (ndata(i) .eq. 0) cycle - scn_strt = scn_end +1 - scn_end = scn_end + ndata(i) - - if (any(flipped_scan2 .eq. remap(i) )) then - do j = scn_strt,scn_end - call flip_e_state(Mux(:,:,j),2) - enddo - endif - if (any(flipped_scan3 .eq. remap(i) )) then - do j = scn_strt,scn_end - call flip_e_state(Mux(:,:,j),3) - enddo - endif - if (any(flipped_scan4 .eq. remap(i))) then - do j =scn_strt,scn_end - call flip_e_state(Mux(:,:,j),4) - enddo - endif - - ii = 4 - do m=1,ndiab - do n = m,ndiab - ii = ii +1 - pat_out(ii,scn_strt:scn_end) = & - Mux(m,n,scn_strt:scn_end) - enddo - enddo - - - - enddo - endif - - deallocate(Mux) - - END SUBROUTINE data_transform - - subroutine Y2mat(Y,Mx) - IMPLICIT NONE - Real(dp), intent(in):: y(:) - Real(dp),intent(out):: Mx(ndiab,ndiab) - ! internal variables - integer(idp):: ii,i,j - - !ntot = ndiab*(ndiab+1)/2 - !if (inp_out/2 .ne. ntot) then - ! write(6,*)"output not equal ntot", inp_out,ntot - ! stop - !endif - ii=1 - do i=1,ndiab - do j=i,ndiab - ! !mx - - mx(i,j)=y(ii) - ! - ii=ii+1 - enddo - enddo - call copy_2_low(mx) - !call copy_2_low(my) - end subroutine - ! ! subroutine which copy the upper matrix to lower - subroutine copy_2_low(m) - implicit none - real(dp), intent(inout) :: m(:,:) - integer(idp) :: i,j - ! copy the lower part of the matrix to the upper part - do i=1,size(m,1) - do j=i,size(m,1) - m(j,i) = m(i,j) - enddo - enddo - end subroutine copy_2_low - - integer function remap(i) - integer(idp),intent(in):: i - - - remap = i + 0 - - end function remap - - - end module data_transf_mod - - diff --git a/src/model/dipole_model.f90 b/src/model/dipole_model.f90 deleted file mode 100644 index 2792908..0000000 --- a/src/model/dipole_model.f90 +++ /dev/null @@ -1,411 +0,0 @@ -module diabmodel - use dip_param - use iso_fortran_env, only:dp => real64,idp => int32 - use nn_params - use nncommons, only: len_in - implicit none - - !integer(idp),parameter:: ndiab=4 - !logical :: debug=.false. - private - public:: diab_x - contains - - - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine diab_x(e,q,nn_out) - real(dp),intent(in) :: q(maxnin) - real(dp),intent(inout) :: nn_out(maxnout) - real(dp),intent(out) :: e(:,:) - integer(idp) id,i,j, ii - real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) - real(dp), dimension(16) :: shift,scal - ! real(dp),dimension(17):: shi,scalee - real(dp), parameter:: tol = 1.0d-9 - !include "shift-scale-70N.incl" - !include "sh-scal-50.incl" - - xs=q(len_in+2) - ys=q(len_in+3) - xb=q(len_in+4) - yb=q(len_in+5) - a=q(len_in+1) - b=q(len_in+6) - - id = 11 - - call init_dip_planar_data() - shift =1.0d0 - scal = 1.0d0 - ii = 5 - do i = 28 ,np ! pes parameter key index - if( abs(p(i)) .gt. tol ) then - p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii)) - ii = ii + 1 - endif - enddo - !write(6,*)" Dipole non-zero", ii-1 - - ss=xs**2+ys**2 ! totaly symmetric term - sb=xb**2+yb**2 - - v3_vec( 1) = xs*(xs**2-3*ys**2) - v3_vec( 2) = xb*(xb**2-3*yb**2) - v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys - v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb - v3_vec( 5) = ys*(3*xs**2-ys**2) - v3_vec( 6) = yb*(3*xb**2-yb**2) - v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys - v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb - - - - e=0.0d0 - - - - !id=key !1 - ! V-term - ! order 1 - e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - id=id+1 !2 - e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - id=id+1 !3 - e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - ! order 2 - id=id+1 !4 - e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb) - id=id+1 !5 - e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +& - p(pst(1,id)+2)*(xs*xb-ys*yb) - e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +& - p(pst(1,id)+2)*(xs*xb-ys*yb) - id=id+1 !6 - e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) + & - p(pst(1,id)+2)*(xs*xb-ys*yb) - ! order 3 - id=id+1 !7 - e(1,1)=e(1,1)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - id=id+1 !8 - e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - id=id+1 !9 - e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb - - ! JAHN TELLER COUPLING W AND Z - ! order 0 - id=id+1 !10 - e(2,2)=e(2,2)+p(pst(1,id)) - e(3,3)=e(3,3)-p(pst(1,id)) - ! order 1 - id=id+1 !11 - e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - e(3,3)=e(3,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb - e(2,3)=e(2,3)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb - ! order 2 - id=id+1 !12 - e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb - e(3,3)=e(3,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb) - e(2,3)=e(2,3)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+ & - p(pst(1,id)+2)*(xs*yb+xb*ys) - ! order 3 - - id=id+1 !13 - do i=1,4 - j=i-1 - e(2,2)=e(2,2)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) - e(3,3)=e(3,3)-(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) - e(2,3)=e(2,3)+(-p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i+4) - enddo - - e(2,2)=e(2,2)+p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb - e(3,3)=e(3,3)-(p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb) - e(2,3)=e(2,3)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb - ! PSEUDO JAHN TELLER - - ! A2 ground state coupled with E - ! ################################################### - ! ################################################### - - ! order 0 - id=id+1 !14 - e(1,2)=e(1,2)+b*p(pst(1,id)) - - ! order 1 - id=id+1 !15 - e(1,2)=e(1,2)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) - e(1,3)=e(1,3)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) - ! order 2 - id=id+1 !16 - e(1,2)=e(1,2)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb) + p(pst(1,id)+3)*(xs**2+ys**2)) - e(1,3)=e(1,3)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& - +p(pst(1,id)+2)*(xs*yb+xb*ys)) - ! order 3 - id =id+1 ! 17 - - do i=1,4 - e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) - e(1,3)=e(1,3)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) - enddo - - - - !! THE COUPLING OF A2 WITH A1 - !#################################################### - !#################################################### - ! order 1 - id=id+1 !18 - e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) - id=id+1 !19 - e(1,4)=e(1,4)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb)) - - - !!! THE COUPLING OF A1 WITH E - !!#################################################### - !#################################################### - ! order 0 - - id=id+1 !20 - e(2,4)=e(2,4)+p(pst(1,id)) - - ! order 1 - id=id+1 !21 - e(2,4)=e(2,4)+p(pst(1,id))*xs + p(pst(1,id)+1)*xb - e(3,4)=e(3,4)+p(pst(1,id))*ys + p(pst(1,id)+1)*yb - - ! order 2 - id=id+1 !22 - e(2,4)=e(2,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb) +p(pst(1,id)+3)*(xs**2+ys**2) - e(3,4)=e(3,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & - -p(pst(1,id)+2)*(xs*yb+xb*ys) - ! order 3 - id=id+1 !23 - do i=1,4 - e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) - e(3,4)=e(3,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) - enddo - - !! End of the model - - e(2,1)=e(1,2) - e(3,1)=e(1,3) - e(3,2)=e(2,3) - e(4,1)=e(1,4) - e(4,2)=e(2,4) - e(4,3)=e(3,4) - end subroutine diab_x - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! THE Y COMPONENT OF DIPOLE - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -!subroutine diab_y(e,q,nn_out) -! real(dp),intent(in) :: q(maxnin) -! real(dp),intent(inout) :: nn_out(maxnout) -! real(dp),intent(out) :: e(:,:) -! integer(idp) id,i,j, ii -! real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) -! real(dp), dimension(23) :: shift,scal -! real(dp), parameter:: tol = 1.0d-9 -! -! xs=q(5) -! ys=q(6) -! xb=q(7) -! yb=q(8) -! a=q(4) -! b=q(9) -! -! id = 11 -! -! call init_dip_planar_data() -! scal = 1.0d-3 -! shift =1.0d0 -! ii = 1 -! do i = 11 ,np ! pes parameter key index -! if( abs(p(i)) .gt. tol ) then -! p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii)) -! ii = ii + 1 -! endif -! enddo -! -! id = 11 -! -! ss=xs**2+ys**2 ! totaly symmetric term -! sb=xb**2+yb**2 -! -! v3_vec( 1) = xs*(xs**2-3*ys**2) -! v3_vec( 2) = xb*(xb**2-3*yb**2) -! v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys -! v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb -! v3_vec( 5) = ys*(3*xs**2-ys**2) -! v3_vec( 6) = yb*(3*xb**2-yb**2) -! v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys -! v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb -! -! -! -! -! -! e=0.0d0 -! ! V-term -! !id=key !1 -! e(1,1)=e(1,1)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! -! id=id+1 !2 -! e(2,2)=e(2,2)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! id=id+1 !3 -! e(4,4)=e(4,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! ! order 2 -! id=id+1 !4 -! e(1,1)=e(1,1)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! id=id+1 !5 -! e(2,2)=e(2,2)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! -! -! e(3,3)=e(3,3)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! id=id+1 !6 -! e(4,4)=e(4,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & -! -p(pst(1,id)+2)*(xs*yb+xb*ys) -! ! order 3 -! id=id+1 !7 -! e(1,1)=e(1,1)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! id=id+1 !8 -! e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! id=id+1 !9 -! e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb -! -! ! V- term + totally symmetric coord a -! -! ! JAHN TELLER COUPLING TERM -! ! order 0 -! id=id+1 !10 -! e(2,3)=e(2,3)+p(pst(1,id)) -! ! order 1 -! -! id=id+1 !11 -! e(2,2)=e(2,2)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb -! e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb -! e(2,3)=e(2,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb -! !id=id+1 !12 -! ! order 2 -! id=id+1 !12 -! e(2,2)=e(2,2)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+p(pst(1,id)+2)*(xs*yb+xb*ys) -! e(3,3)=e(3,3)-p(pst(1,id))*2*xs*ys-p(pst(1,id)+1)*2*xb*yb-p(pst(1,id)+2)*(xs*yb+xb*ys) -! e(2,3)=e(2,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)) & -! -p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb -! ! order 3 -! id=id+1 !13 -! do i=1,4 -! j=i-1 -! e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) -! e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) -! e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) -! enddo -! e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb -! e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb -! e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb -! ! PSEUDO JAHN TELLER -! ! ORDER 0 -! ! THE COUPLING OF A2 GROUND STATE WITH E -! ! ################################################### -! ! ################################################### -! ! order 0 -! id=id+1 !14 -! e(1,3)=e(1,3)-b*(p(pst(1,id))) -! ! order 1 -! id=id+1 !15 -! e(1,2)=e(1,2)-b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) -! e(1,3)=e(1,3)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) -! -! ! order 2 -! id=id+1 !16 -! e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& -! +p(pst(1,id)+2)*(xs*yb+xb*ys)) -! e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& -! +p(pst(1,id)+2)*(xs*xb-ys*yb) - p(pst(1,id)+3)*(xs**2+ys**2)) -! -! ! order 3 -! id = id+1 ! 17 -! do i=1,4 -! e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) -! e(1,3)=e(1,3)-b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) -! enddo -! -! -! ! THE COUPLING OF A2 WITH A1 -! !#################################################### -! !#################################################### -! ! order 1 -! id=id+1 !17 -! e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) -! ! order 2 -! id=id+1 !18 -! e(1,4)=e(1,4)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& -! +p(pst(1,id)+2)*(xs*yb+xb*ys)) -! -! -! ! THE COUPLING OF A1 WITH E -! !#################################################### -! !#################################################### -! ! order 0 -! id=id+1 !19 -! e(3,4)=e(3,4)-p(pst(1,id)) -! ! order 1 -! id=id+1 !20 -! e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb -! e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb -! ! order 2 -! id=id+1 !21 -! e(2,4)=e(2,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb) & -! +p(pst(1,id)+2)*(xs*yb+xb*ys) -! e(3,4)=e(3,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & -! +p(pst(1,id)+2)*(xs*xb-ys*yb) - p(pst(1,id)+3)*(xs**2+ys**2) -! -! id =id+1 ! 23 -! ! order 3 -! do i=1,4 -! e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) -! e(3,4)=e(3,4)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) -! enddo -! -! ! end of the model -! e(2,1)=e(1,2) -! e(3,1)=e(1,3) -! e(3,2)=e(2,3) -! e(4,1)=e(1,4) -! e(4,2)=e(2,4) -! e(4,3)=e(3,4) -! ! change Ex and Ey -!end subroutine diab_y - -! potential nh3 -! parth - - -!subroutine copy_2_lower_triangle(mat) -! real(dp), intent(inout) :: mat(:, :) -! integer :: m, n -! write lower triangle of matrix symmetrical -! do n = 2, size(mat, 1) -! do m = 1, n - 1 -! mat(n, m) = mat(m, n) -! end do -! end do -! end subroutine copy_2_lower_triangle - -end module diabmodel diff --git a/src/model/fit_genetic.f90 b/src/model/fit_genetic.f90 deleted file mode 100644 index b1c0bc0..0000000 --- a/src/model/fit_genetic.f90 +++ /dev/null @@ -1,914 +0,0 @@ - Module dip_param - IMPLICIT NONE - Integer,parameter :: np=101 - Double precision :: p(101) - integer :: pst(2,400) - ! Non zero parameters are 17 - !$omp threadprivate(p,pst) - contains - - SUBROUTINE init_dip_planar_data() - implicit none - p( 1)= -0.561433980D+02 - p( 2)= 0.000000000D+00 - p( 3)= 0.000000000D+00 - p( 4)= -0.558830090D+02 - p( 5)= 0.000000000D+00 - p( 6)= 0.000000000D+00 - p( 7)= -0.557220030D+02 - p( 8)= 0.000000000D+00 - p( 9)= 0.000000000D+00 - p( 10)= 0.000000000D+00 - p( 11)= 0.334452655D-02 - p( 12)= 0.000000000D+00 - p( 13)= 0.000000000D+00 - p( 14)= 0.000000000D+00 - p( 15)= 0.000000000D+00 - p( 16)= 0.000000000D+00 - p( 17)= 0.000000000D+00 - p( 18)= 0.000000000D+00 - p( 19)= 0.000000000D+00 - p( 20)= 0.000000000D+00 - p( 21)= -0.532590808D-03 - p( 22)= 0.000000000D+00 - p( 23)= 0.000000000D+00 - p( 24)= 0.000000000D+00 - p( 25)= 0.000000000D+00 - p( 26)= 0.000000000D+00 - p( 27)= 0.000000000D+00 - p( 28)= -0.226116045D-01 - p( 29)= 0.000000000D+00 - p( 30)= -0.363599361D-01 - p( 31)= 0.000000000D+00 - p( 32)= 0.962100323D-01 - p( 33)= 0.000000000D+00 - p( 34)= 0.111336127D-03 - p( 35)= 0.000000000D+00 - p( 36)= 0.000000000D+00 - p( 37)= -0.433685770D-04 - p( 38)= 0.000000000D+00 - p( 39)= 0.000000000D+00 - p( 40)= -0.548870572D-02 - p( 41)= 0.000000000D+00 - p( 42)= 0.000000000D+00 - p( 43)= 0.000000000D+00 - p( 44)= 0.000000000D+00 - p( 45)= 0.000000000D+00 - p( 46)= 0.000000000D+00 - p( 47)= 0.000000000D+00 - p( 48)= 0.000000000D+00 - p( 49)= 0.217598103D+00 - p( 50)= -0.113909871D-01 - p( 51)= 0.000000000D+00 - p( 52)= -0.123184224D-03 - p( 53)= 0.000000000D+00 - p( 54)= 0.000000000D+00 - p( 55)= 0.000000000D+00 - p( 56)= 0.000000000D+00 - p( 57)= 0.000000000D+00 - p( 58)= 0.000000000D+00 - p( 59)= 0.000000000D+00 - p( 60)= 0.000000000D+00 - p( 61)= 0.000000000D+00 - p( 62)= 0.000000000D+00 - p( 63)= 0.000000000D+00 - p( 64)= 0.000000000D+00 - p( 65)= 0.000000000D+00 - p( 66)= 0.000000000D+00 - p( 67)= 0.000000000D+00 - p( 68)= 0.000000000D+00 - p( 69)= 0.000000000D+00 - p( 70)= 0.000000000D+00 - p( 71)= 0.000000000D+00 - p( 72)= 0.000000000D+00 - p( 73)= 0.000000000D+00 - p( 74)= 0.000000000D+00 - p( 75)= 0.000000000D+00 - p( 76)= 0.000000000D+00 - p( 77)= 0.000000000D+00 - p( 78)= 0.000000000D+00 - p( 79)= 0.000000000D+00 - p( 80)= 0.000000000D+00 - p( 81)= 0.000000000D+00 - p( 82)= 0.000000000D+00 - p( 83)= 0.000000000D+00 - p( 84)= 0.000000000D+00 - p( 85)= 0.000000000D+00 - p( 86)= 0.000000000D+00 - p( 87)= -0.567206500D-02 - p( 88)= -0.613490432D-02 - p( 89)= 0.000000000D+00 - p( 90)= 0.303022045D-03 - p( 91)= 0.000000000D+00 - p( 92)= 0.000000000D+00 - p( 93)= 0.000000000D+00 - p( 94)= 0.000000000D+00 - p( 95)= 0.000000000D+00 - p( 96)= 0.000000000D+00 - p( 97)= 0.000000000D+00 - p( 98)= 0.000000000D+00 - p( 99)= 0.000000000D+00 - p(100)= 0.000000000D+00 - p(101)= 0.000000000D+00 - pst(1, 1)= 1 - pst(2, 1)= 3 - pst(1, 2)= 4 - pst(2, 2)= 3 - pst(1, 3)= 7 - pst(2, 3)= 4 - pst(1, 4)= 11 - pst(2, 4)= 2 - pst(1, 5)= 13 - pst(2, 5)= 3 - pst(1, 6)= 16 - pst(2, 6)= 2 - pst(1, 7)= 18 - pst(2, 7)= 3 - pst(1, 8)= 21 - pst(2, 8)= 3 - pst(1, 9)= 24 - pst(2, 9)= 3 - pst(1, 10)= 27 - pst(2, 10)= 1 - pst(1, 11)= 28 - pst(2, 11)= 2 - pst(1, 12)= 30 - pst(2, 12)= 2 - pst(1, 13)= 32 - pst(2, 13)= 2 - pst(1, 14)= 34 - pst(2, 14)= 3 - pst(1, 15)= 37 - pst(2, 15)= 3 - pst(1, 16)= 40 - pst(2, 16)= 3 - pst(1, 17)= 43 - pst(2, 17)= 2 - pst(1, 18)= 45 - pst(2, 18)= 2 - pst(1, 19)= 47 - pst(2, 19)= 2 - pst(1, 20)= 49 - pst(2, 20)= 1 - pst(1, 21)= 50 - pst(2, 21)= 2 - pst(1, 22)= 52 - pst(2, 22)= 5 - pst(1, 23)= 57 - pst(2, 23)= 10 - pst(1, 24)= 67 - pst(2, 24)= 1 - pst(1, 25)= 68 - pst(2, 25)= 2 - pst(1, 26)= 70 - pst(2, 26)= 4 - pst(1, 27)= 74 - pst(2, 27)= 8 - pst(1, 28)= 82 - pst(2, 28)= 2 - pst(1, 29)= 84 - pst(2, 29)= 3 - pst(1, 30)= 87 - pst(2, 30)= 1 - pst(1, 31)= 88 - pst(2, 31)= 2 - pst(1, 32)= 90 - pst(2, 32)= 4 - pst(1, 33)= 94 - pst(2, 33)= 8 - pst(1, 34)= 0 - pst(2, 34)= 0 - pst(1, 35)= 0 - pst(2, 35)= 0 - pst(1, 36)= 0 - pst(2, 36)= 0 - pst(1, 37)= 0 - pst(2, 37)= 0 - pst(1, 38)= 0 - pst(2, 38)= 0 - pst(1, 39)= 0 - pst(2, 39)= 0 - pst(1, 40)= 0 - pst(2, 40)= 0 - pst(1, 41)= 0 - pst(2, 41)= 0 - pst(1, 42)= 0 - pst(2, 42)= 0 - pst(1, 43)= 0 - pst(2, 43)= 0 - pst(1, 44)= 0 - pst(2, 44)= 0 - pst(1, 45)= 0 - pst(2, 45)= 0 - pst(1, 46)= 0 - pst(2, 46)= 0 - pst(1, 47)= 0 - pst(2, 47)= 0 - pst(1, 48)= 0 - pst(2, 48)= 0 - pst(1, 49)= 0 - pst(2, 49)= 0 - pst(1, 50)= 0 - pst(2, 50)= 0 - pst(1, 51)= 0 - pst(2, 51)= 0 - pst(1, 52)= 0 - pst(2, 52)= 0 - pst(1, 53)= 0 - pst(2, 53)= 0 - pst(1, 54)= 0 - pst(2, 54)= 0 - pst(1, 55)= 0 - pst(2, 55)= 0 - pst(1, 56)= 0 - pst(2, 56)= 0 - pst(1, 57)= 0 - pst(2, 57)= 0 - pst(1, 58)= 0 - pst(2, 58)= 0 - pst(1, 59)= 0 - pst(2, 59)= 0 - pst(1, 60)= 0 - pst(2, 60)= 0 - pst(1, 61)= 0 - pst(2, 61)= 0 - pst(1, 62)= 0 - pst(2, 62)= 0 - pst(1, 63)= 0 - pst(2, 63)= 0 - pst(1, 64)= 0 - pst(2, 64)= 0 - pst(1, 65)= 0 - pst(2, 65)= 0 - pst(1, 66)= 0 - pst(2, 66)= 0 - pst(1, 67)= 0 - pst(2, 67)= 0 - pst(1, 68)= 0 - pst(2, 68)= 0 - pst(1, 69)= 0 - pst(2, 69)= 0 - pst(1, 70)= 0 - pst(2, 70)= 0 - pst(1, 71)= 0 - pst(2, 71)= 0 - pst(1, 72)= 0 - pst(2, 72)= 0 - pst(1, 73)= 0 - pst(2, 73)= 0 - pst(1, 74)= 0 - pst(2, 74)= 0 - pst(1, 75)= 0 - pst(2, 75)= 0 - pst(1, 76)= 0 - pst(2, 76)= 0 - pst(1, 77)= 0 - pst(2, 77)= 0 - pst(1, 78)= 0 - pst(2, 78)= 0 - pst(1, 79)= 0 - pst(2, 79)= 0 - pst(1, 80)= 0 - pst(2, 80)= 0 - pst(1, 81)= 0 - pst(2, 81)= 0 - pst(1, 82)= 0 - pst(2, 82)= 0 - pst(1, 83)= 0 - pst(2, 83)= 0 - pst(1, 84)= 0 - pst(2, 84)= 0 - pst(1, 85)= 0 - pst(2, 85)= 0 - pst(1, 86)= 0 - pst(2, 86)= 0 - pst(1, 87)= 0 - pst(2, 87)= 0 - pst(1, 88)= 0 - pst(2, 88)= 0 - pst(1, 89)= 0 - pst(2, 89)= 0 - pst(1, 90)= 0 - pst(2, 90)= 0 - pst(1, 91)= 0 - pst(2, 91)= 0 - pst(1, 92)= 0 - pst(2, 92)= 0 - pst(1, 93)= 0 - pst(2, 93)= 0 - pst(1, 94)= 0 - pst(2, 94)= 0 - pst(1, 95)= 0 - pst(2, 95)= 0 - pst(1, 96)= 0 - pst(2, 96)= 0 - pst(1, 97)= 0 - pst(2, 97)= 0 - pst(1, 98)= 0 - pst(2, 98)= 0 - pst(1, 99)= 0 - pst(2, 99)= 0 - pst(1,100)= 0 - pst(2,100)= 0 - pst(1,101)= 0 - pst(2,101)= 0 - pst(1,102)= 0 - pst(2,102)= 0 - pst(1,103)= 0 - pst(2,103)= 0 - pst(1,104)= 0 - pst(2,104)= 0 - pst(1,105)= 0 - pst(2,105)= 0 - pst(1,106)= 0 - pst(2,106)= 0 - pst(1,107)= 0 - pst(2,107)= 0 - pst(1,108)= 0 - pst(2,108)= 0 - pst(1,109)= 0 - pst(2,109)= 0 - pst(1,110)= 0 - pst(2,110)= 0 - pst(1,111)= 0 - pst(2,111)= 0 - pst(1,112)= 0 - pst(2,112)= 0 - pst(1,113)= 0 - pst(2,113)= 0 - pst(1,114)= 0 - pst(2,114)= 0 - pst(1,115)= 0 - pst(2,115)= 0 - pst(1,116)= 0 - pst(2,116)= 0 - pst(1,117)= 0 - pst(2,117)= 0 - pst(1,118)= 0 - pst(2,118)= 0 - pst(1,119)= 0 - pst(2,119)= 0 - pst(1,120)= 0 - pst(2,120)= 0 - pst(1,121)= 0 - pst(2,121)= 0 - pst(1,122)= 0 - pst(2,122)= 0 - pst(1,123)= 0 - pst(2,123)= 0 - pst(1,124)= 0 - pst(2,124)= 0 - pst(1,125)= 0 - pst(2,125)= 0 - pst(1,126)= 0 - pst(2,126)= 0 - pst(1,127)= 0 - pst(2,127)= 0 - pst(1,128)= 0 - pst(2,128)= 0 - pst(1,129)= 0 - pst(2,129)= 0 - pst(1,130)= 0 - pst(2,130)= 0 - pst(1,131)= 0 - pst(2,131)= 0 - pst(1,132)= 0 - pst(2,132)= 0 - pst(1,133)= 0 - pst(2,133)= 0 - pst(1,134)= 0 - pst(2,134)= 0 - pst(1,135)= 0 - pst(2,135)= 0 - pst(1,136)= 0 - pst(2,136)= 0 - pst(1,137)= 0 - pst(2,137)= 0 - pst(1,138)= 0 - pst(2,138)= 0 - pst(1,139)= 0 - pst(2,139)= 0 - pst(1,140)= 0 - pst(2,140)= 0 - pst(1,141)= 0 - pst(2,141)= 0 - pst(1,142)= 0 - pst(2,142)= 0 - pst(1,143)= 0 - pst(2,143)= 0 - pst(1,144)= 0 - pst(2,144)= 0 - pst(1,145)= 0 - pst(2,145)= 0 - pst(1,146)= 0 - pst(2,146)= 0 - pst(1,147)= 0 - pst(2,147)= 0 - pst(1,148)= 0 - pst(2,148)= 0 - pst(1,149)= 0 - pst(2,149)= 0 - pst(1,150)= 0 - pst(2,150)= 0 - pst(1,151)= 0 - pst(2,151)= 0 - pst(1,152)= 0 - pst(2,152)= 0 - pst(1,153)= 0 - pst(2,153)= 0 - pst(1,154)= 0 - pst(2,154)= 0 - pst(1,155)= 0 - pst(2,155)= 0 - pst(1,156)= 0 - pst(2,156)= 0 - pst(1,157)= 0 - pst(2,157)= 0 - pst(1,158)= 0 - pst(2,158)= 0 - pst(1,159)= 0 - pst(2,159)= 0 - pst(1,160)= 0 - pst(2,160)= 0 - pst(1,161)= 0 - pst(2,161)= 0 - pst(1,162)= 0 - pst(2,162)= 0 - pst(1,163)= 0 - pst(2,163)= 0 - pst(1,164)= 0 - pst(2,164)= 0 - pst(1,165)= 0 - pst(2,165)= 0 - pst(1,166)= 0 - pst(2,166)= 0 - pst(1,167)= 0 - pst(2,167)= 0 - pst(1,168)= 0 - pst(2,168)= 0 - pst(1,169)= 0 - pst(2,169)= 0 - pst(1,170)= 0 - pst(2,170)= 0 - pst(1,171)= 0 - pst(2,171)= 0 - pst(1,172)= 0 - pst(2,172)= 0 - pst(1,173)= 0 - pst(2,173)= 0 - pst(1,174)= 0 - pst(2,174)= 0 - pst(1,175)= 0 - pst(2,175)= 0 - pst(1,176)= 0 - pst(2,176)= 0 - pst(1,177)= 0 - pst(2,177)= 0 - pst(1,178)= 0 - pst(2,178)= 0 - pst(1,179)= 0 - pst(2,179)= 0 - pst(1,180)= 0 - pst(2,180)= 0 - pst(1,181)= 0 - pst(2,181)= 0 - pst(1,182)= 0 - pst(2,182)= 0 - pst(1,183)= 0 - pst(2,183)= 0 - pst(1,184)= 0 - pst(2,184)= 0 - pst(1,185)= 0 - pst(2,185)= 0 - pst(1,186)= 0 - pst(2,186)= 0 - pst(1,187)= 0 - pst(2,187)= 0 - pst(1,188)= 0 - pst(2,188)= 0 - pst(1,189)= 0 - pst(2,189)= 0 - pst(1,190)= 0 - pst(2,190)= 0 - pst(1,191)= 0 - pst(2,191)= 0 - pst(1,192)= 0 - pst(2,192)= 0 - pst(1,193)= 0 - pst(2,193)= 0 - pst(1,194)= 0 - pst(2,194)= 0 - pst(1,195)= 0 - pst(2,195)= 0 - pst(1,196)= 0 - pst(2,196)= 0 - pst(1,197)= 0 - pst(2,197)= 0 - pst(1,198)= 0 - pst(2,198)= 0 - pst(1,199)= 0 - pst(2,199)= 0 - pst(1,200)= 0 - pst(2,200)= 0 - pst(1,201)= 0 - pst(2,201)= 0 - pst(1,202)= 0 - pst(2,202)= 0 - pst(1,203)= 0 - pst(2,203)= 0 - pst(1,204)= 0 - pst(2,204)= 0 - pst(1,205)= 0 - pst(2,205)= 0 - pst(1,206)= 0 - pst(2,206)= 0 - pst(1,207)= 0 - pst(2,207)= 0 - pst(1,208)= 0 - pst(2,208)= 0 - pst(1,209)= 0 - pst(2,209)= 0 - pst(1,210)= 0 - pst(2,210)= 0 - pst(1,211)= 0 - pst(2,211)= 0 - pst(1,212)= 0 - pst(2,212)= 0 - pst(1,213)= 0 - pst(2,213)= 0 - pst(1,214)= 0 - pst(2,214)= 0 - pst(1,215)= 0 - pst(2,215)= 0 - pst(1,216)= 0 - pst(2,216)= 0 - pst(1,217)= 0 - pst(2,217)= 0 - pst(1,218)= 0 - pst(2,218)= 0 - pst(1,219)= 0 - pst(2,219)= 0 - pst(1,220)= 0 - pst(2,220)= 0 - pst(1,221)= 0 - pst(2,221)= 0 - pst(1,222)= 0 - pst(2,222)= 0 - pst(1,223)= 0 - pst(2,223)= 0 - pst(1,224)= 0 - pst(2,224)= 0 - pst(1,225)= 0 - pst(2,225)= 0 - pst(1,226)= 0 - pst(2,226)= 0 - pst(1,227)= 0 - pst(2,227)= 0 - pst(1,228)= 0 - pst(2,228)= 0 - pst(1,229)= 0 - pst(2,229)= 0 - pst(1,230)= 0 - pst(2,230)= 0 - pst(1,231)= 0 - pst(2,231)= 0 - pst(1,232)= 0 - pst(2,232)= 0 - pst(1,233)= 0 - pst(2,233)= 0 - pst(1,234)= 0 - pst(2,234)= 0 - pst(1,235)= 0 - pst(2,235)= 0 - pst(1,236)= 0 - pst(2,236)= 0 - pst(1,237)= 0 - pst(2,237)= 0 - pst(1,238)= 0 - pst(2,238)= 0 - pst(1,239)= 0 - pst(2,239)= 0 - pst(1,240)= 0 - pst(2,240)= 0 - pst(1,241)= 0 - pst(2,241)= 0 - pst(1,242)= 0 - pst(2,242)= 0 - pst(1,243)= 0 - pst(2,243)= 0 - pst(1,244)= 0 - pst(2,244)= 0 - pst(1,245)= 0 - pst(2,245)= 0 - pst(1,246)= 0 - pst(2,246)= 0 - pst(1,247)= 0 - pst(2,247)= 0 - pst(1,248)= 0 - pst(2,248)= 0 - pst(1,249)= 0 - pst(2,249)= 0 - pst(1,250)= 0 - pst(2,250)= 0 - pst(1,251)= 0 - pst(2,251)= 0 - pst(1,252)= 0 - pst(2,252)= 0 - pst(1,253)= 0 - pst(2,253)= 0 - pst(1,254)= 0 - pst(2,254)= 0 - pst(1,255)= 0 - pst(2,255)= 0 - pst(1,256)= 0 - pst(2,256)= 0 - pst(1,257)= 0 - pst(2,257)= 0 - pst(1,258)= 0 - pst(2,258)= 0 - pst(1,259)= 0 - pst(2,259)= 0 - pst(1,260)= 0 - pst(2,260)= 0 - pst(1,261)= 0 - pst(2,261)= 0 - pst(1,262)= 0 - pst(2,262)= 0 - pst(1,263)= 0 - pst(2,263)= 0 - pst(1,264)= 0 - pst(2,264)= 0 - pst(1,265)= 0 - pst(2,265)= 0 - pst(1,266)= 0 - pst(2,266)= 0 - pst(1,267)= 0 - pst(2,267)= 0 - pst(1,268)= 0 - pst(2,268)= 0 - pst(1,269)= 0 - pst(2,269)= 0 - pst(1,270)= 0 - pst(2,270)= 0 - pst(1,271)= 0 - pst(2,271)= 0 - pst(1,272)= 0 - pst(2,272)= 0 - pst(1,273)= 0 - pst(2,273)= 0 - pst(1,274)= 0 - pst(2,274)= 0 - pst(1,275)= 0 - pst(2,275)= 0 - pst(1,276)= 0 - pst(2,276)= 0 - pst(1,277)= 0 - pst(2,277)= 0 - pst(1,278)= 0 - pst(2,278)= 0 - pst(1,279)= 0 - pst(2,279)= 0 - pst(1,280)= 0 - pst(2,280)= 0 - pst(1,281)= 0 - pst(2,281)= 0 - pst(1,282)= 0 - pst(2,282)= 0 - pst(1,283)= 0 - pst(2,283)= 0 - pst(1,284)= 0 - pst(2,284)= 0 - pst(1,285)= 0 - pst(2,285)= 0 - pst(1,286)= 0 - pst(2,286)= 0 - pst(1,287)= 0 - pst(2,287)= 0 - pst(1,288)= 0 - pst(2,288)= 0 - pst(1,289)= 0 - pst(2,289)= 0 - pst(1,290)= 0 - pst(2,290)= 0 - pst(1,291)= 0 - pst(2,291)= 0 - pst(1,292)= 0 - pst(2,292)= 0 - pst(1,293)= 0 - pst(2,293)= 0 - pst(1,294)= 0 - pst(2,294)= 0 - pst(1,295)= 0 - pst(2,295)= 0 - pst(1,296)= 0 - pst(2,296)= 0 - pst(1,297)= 0 - pst(2,297)= 0 - pst(1,298)= 0 - pst(2,298)= 0 - pst(1,299)= 0 - pst(2,299)= 0 - pst(1,300)= 0 - pst(2,300)= 0 - pst(1,301)= 0 - pst(2,301)= 0 - pst(1,302)= 0 - pst(2,302)= 0 - pst(1,303)= 0 - pst(2,303)= 0 - pst(1,304)= 0 - pst(2,304)= 0 - pst(1,305)= 0 - pst(2,305)= 0 - pst(1,306)= 0 - pst(2,306)= 0 - pst(1,307)= 0 - pst(2,307)= 0 - pst(1,308)= 0 - pst(2,308)= 0 - pst(1,309)= 0 - pst(2,309)= 0 - pst(1,310)= 0 - pst(2,310)= 0 - pst(1,311)= 0 - pst(2,311)= 0 - pst(1,312)= 0 - pst(2,312)= 0 - pst(1,313)= 0 - pst(2,313)= 0 - pst(1,314)= 0 - pst(2,314)= 0 - pst(1,315)= 0 - pst(2,315)= 0 - pst(1,316)= 0 - pst(2,316)= 0 - pst(1,317)= 0 - pst(2,317)= 0 - pst(1,318)= 0 - pst(2,318)= 0 - pst(1,319)= 0 - pst(2,319)= 0 - pst(1,320)= 0 - pst(2,320)= 0 - pst(1,321)= 0 - pst(2,321)= 0 - pst(1,322)= 0 - pst(2,322)= 0 - pst(1,323)= 0 - pst(2,323)= 0 - pst(1,324)= 0 - pst(2,324)= 0 - pst(1,325)= 0 - pst(2,325)= 0 - pst(1,326)= 0 - pst(2,326)= 0 - pst(1,327)= 0 - pst(2,327)= 0 - pst(1,328)= 0 - pst(2,328)= 0 - pst(1,329)= 0 - pst(2,329)= 0 - pst(1,330)= 0 - pst(2,330)= 0 - pst(1,331)= 0 - pst(2,331)= 0 - pst(1,332)= 0 - pst(2,332)= 0 - pst(1,333)= 0 - pst(2,333)= 0 - pst(1,334)= 0 - pst(2,334)= 0 - pst(1,335)= 0 - pst(2,335)= 0 - pst(1,336)= 0 - pst(2,336)= 0 - pst(1,337)= 0 - pst(2,337)= 0 - pst(1,338)= 0 - pst(2,338)= 0 - pst(1,339)= 0 - pst(2,339)= 0 - pst(1,340)= 0 - pst(2,340)= 0 - pst(1,341)= 0 - pst(2,341)= 0 - pst(1,342)= 0 - pst(2,342)= 0 - pst(1,343)= 0 - pst(2,343)= 0 - pst(1,344)= 0 - pst(2,344)= 0 - pst(1,345)= 0 - pst(2,345)= 0 - pst(1,346)= 0 - pst(2,346)= 0 - pst(1,347)= 0 - pst(2,347)= 0 - pst(1,348)= 0 - pst(2,348)= 0 - pst(1,349)= 0 - pst(2,349)= 0 - pst(1,350)= 0 - pst(2,350)= 0 - pst(1,351)= 0 - pst(2,351)= 0 - pst(1,352)= 0 - pst(2,352)= 0 - pst(1,353)= 0 - pst(2,353)= 0 - pst(1,354)= 0 - pst(2,354)= 0 - pst(1,355)= 0 - pst(2,355)= 0 - pst(1,356)= 0 - pst(2,356)= 0 - pst(1,357)= 0 - pst(2,357)= 0 - pst(1,358)= 0 - pst(2,358)= 0 - pst(1,359)= 0 - pst(2,359)= 0 - pst(1,360)= 0 - pst(2,360)= 0 - pst(1,361)= 0 - pst(2,361)= 0 - pst(1,362)= 0 - pst(2,362)= 0 - pst(1,363)= 0 - pst(2,363)= 0 - pst(1,364)= 0 - pst(2,364)= 0 - pst(1,365)= 0 - pst(2,365)= 0 - pst(1,366)= 0 - pst(2,366)= 0 - pst(1,367)= 0 - pst(2,367)= 0 - pst(1,368)= 0 - pst(2,368)= 0 - pst(1,369)= 0 - pst(2,369)= 0 - pst(1,370)= 0 - pst(2,370)= 0 - pst(1,371)= 0 - pst(2,371)= 0 - pst(1,372)= 0 - pst(2,372)= 0 - pst(1,373)= 0 - pst(2,373)= 0 - pst(1,374)= 0 - pst(2,374)= 0 - pst(1,375)= 0 - pst(2,375)= 0 - pst(1,376)= 0 - pst(2,376)= 0 - pst(1,377)= 0 - pst(2,377)= 0 - pst(1,378)= 0 - pst(2,378)= 0 - pst(1,379)= 0 - pst(2,379)= 0 - pst(1,380)= 0 - pst(2,380)= 0 - pst(1,381)= 0 - pst(2,381)= 0 - pst(1,382)= 0 - pst(2,382)= 0 - pst(1,383)= 0 - pst(2,383)= 0 - pst(1,384)= 0 - pst(2,384)= 0 - pst(1,385)= 0 - pst(2,385)= 0 - pst(1,386)= 0 - pst(2,386)= 0 - pst(1,387)= 0 - pst(2,387)= 0 - pst(1,388)= 0 - pst(2,388)= 0 - pst(1,389)= 0 - pst(2,389)= 0 - pst(1,390)= 0 - pst(2,390)= 0 - pst(1,391)= 0 - pst(2,391)= 0 - pst(1,392)= 0 - pst(2,392)= 0 - pst(1,393)= 0 - pst(2,393)= 0 - pst(1,394)= 0 - pst(2,394)= 0 - pst(1,395)= 0 - pst(2,395)= 0 - pst(1,396)= 0 - pst(2,396)= 0 - pst(1,397)= 0 - pst(2,397)= 0 - pst(1,398)= 0 - pst(2,398)= 0 - pst(1,399)= 0 - pst(2,399)= 0 - pst(1,400)= 0 - pst(2,400)= 0 -End SUBROUTINE init_dip_planar_data -End Module dip_param diff --git a/src/model/fit_param_dip_pes.f90 b/src/model/fit_param_dip_pes.f90 deleted file mode 100644 index 372283a..0000000 --- a/src/model/fit_param_dip_pes.f90 +++ /dev/null @@ -1,914 +0,0 @@ - Module dip_param - IMPLICIT NONE - Integer,parameter :: np=101 - Double precision :: p(101) - integer :: pst(2,400) - ! Non zero parameters are 16 - !$omp threadprivate(p,pst) - contains - - SUBROUTINE init_dip_planar_data() - implicit none - p( 1)= -0.561433980D+02 - p( 2)= 0.000000000D+00 - p( 3)= 0.000000000D+00 - p( 4)= -0.558830090D+02 - p( 5)= 0.000000000D+00 - p( 6)= 0.000000000D+00 - p( 7)= -0.557220030D+02 - p( 8)= 0.000000000D+00 - p( 9)= 0.000000000D+00 - p( 10)= 0.000000000D+00 - p( 11)= 0.333624014D-02 - p( 12)= 0.000000000D+00 - p( 13)= 0.000000000D+00 - p( 14)= 0.000000000D+00 - p( 15)= 0.000000000D+00 - p( 16)= 0.000000000D+00 - p( 17)= 0.000000000D+00 - p( 18)= 0.000000000D+00 - p( 19)= 0.000000000D+00 - p( 20)= 0.000000000D+00 - p( 21)= -0.532590808D-03 - p( 22)= 0.000000000D+00 - p( 23)= 0.000000000D+00 - p( 24)= 0.000000000D+00 - p( 25)= 0.000000000D+00 - p( 26)= 0.000000000D+00 - p( 27)= 0.000000000D+00 - p( 28)= -0.226073502D-01 - p( 29)= 0.000000000D+00 - p( 30)= -0.363395355D-01 - p( 31)= 0.000000000D+00 - p( 32)= 0.939930347D-01 - p( 33)= 0.000000000D+00 - p( 34)= 0.111364489D-03 - p( 35)= 0.000000000D+00 - p( 36)= 0.000000000D+00 - p( 37)= 0.000000000D+00 - p( 38)= 0.000000000D+00 - p( 39)= 0.000000000D+00 - p( 40)= -0.545384229D-02 - p( 41)= 0.000000000D+00 - p( 42)= 0.000000000D+00 - p( 43)= 0.000000000D+00 - p( 44)= 0.000000000D+00 - p( 45)= 0.000000000D+00 - p( 46)= 0.000000000D+00 - p( 47)= 0.000000000D+00 - p( 48)= 0.000000000D+00 - p( 49)= 0.217598103D+00 - p( 50)= -0.113772250D-01 - p( 51)= 0.000000000D+00 - p( 52)= -0.122160662D-03 - p( 53)= 0.000000000D+00 - p( 54)= 0.000000000D+00 - p( 55)= 0.000000000D+00 - p( 56)= 0.000000000D+00 - p( 57)= 0.000000000D+00 - p( 58)= 0.000000000D+00 - p( 59)= 0.000000000D+00 - p( 60)= 0.000000000D+00 - p( 61)= 0.000000000D+00 - p( 62)= 0.000000000D+00 - p( 63)= 0.000000000D+00 - p( 64)= 0.000000000D+00 - p( 65)= 0.000000000D+00 - p( 66)= 0.000000000D+00 - p( 67)= 0.000000000D+00 - p( 68)= 0.000000000D+00 - p( 69)= 0.000000000D+00 - p( 70)= 0.000000000D+00 - p( 71)= 0.000000000D+00 - p( 72)= 0.000000000D+00 - p( 73)= 0.000000000D+00 - p( 74)= 0.000000000D+00 - p( 75)= 0.000000000D+00 - p( 76)= 0.000000000D+00 - p( 77)= 0.000000000D+00 - p( 78)= 0.000000000D+00 - p( 79)= 0.000000000D+00 - p( 80)= 0.000000000D+00 - p( 81)= 0.000000000D+00 - p( 82)= 0.000000000D+00 - p( 83)= 0.000000000D+00 - p( 84)= 0.000000000D+00 - p( 85)= 0.000000000D+00 - p( 86)= 0.000000000D+00 - p( 87)= -0.567206500D-02 - p( 88)= -0.565873672D-02 - p( 89)= 0.000000000D+00 - p( 90)= 0.336524505D-03 - p( 91)= 0.000000000D+00 - p( 92)= 0.000000000D+00 - p( 93)= 0.000000000D+00 - p( 94)= 0.000000000D+00 - p( 95)= 0.000000000D+00 - p( 96)= 0.000000000D+00 - p( 97)= 0.000000000D+00 - p( 98)= 0.000000000D+00 - p( 99)= 0.000000000D+00 - p(100)= 0.000000000D+00 - p(101)= 0.000000000D+00 - pst(1, 1)= 1 - pst(2, 1)= 3 - pst(1, 2)= 4 - pst(2, 2)= 3 - pst(1, 3)= 7 - pst(2, 3)= 4 - pst(1, 4)= 11 - pst(2, 4)= 2 - pst(1, 5)= 13 - pst(2, 5)= 3 - pst(1, 6)= 16 - pst(2, 6)= 2 - pst(1, 7)= 18 - pst(2, 7)= 3 - pst(1, 8)= 21 - pst(2, 8)= 3 - pst(1, 9)= 24 - pst(2, 9)= 3 - pst(1, 10)= 27 - pst(2, 10)= 1 - pst(1, 11)= 28 - pst(2, 11)= 2 - pst(1, 12)= 30 - pst(2, 12)= 2 - pst(1, 13)= 32 - pst(2, 13)= 2 - pst(1, 14)= 34 - pst(2, 14)= 3 - pst(1, 15)= 37 - pst(2, 15)= 3 - pst(1, 16)= 40 - pst(2, 16)= 3 - pst(1, 17)= 43 - pst(2, 17)= 2 - pst(1, 18)= 45 - pst(2, 18)= 2 - pst(1, 19)= 47 - pst(2, 19)= 2 - pst(1, 20)= 49 - pst(2, 20)= 1 - pst(1, 21)= 50 - pst(2, 21)= 2 - pst(1, 22)= 52 - pst(2, 22)= 5 - pst(1, 23)= 57 - pst(2, 23)= 10 - pst(1, 24)= 67 - pst(2, 24)= 1 - pst(1, 25)= 68 - pst(2, 25)= 2 - pst(1, 26)= 70 - pst(2, 26)= 4 - pst(1, 27)= 74 - pst(2, 27)= 8 - pst(1, 28)= 82 - pst(2, 28)= 2 - pst(1, 29)= 84 - pst(2, 29)= 3 - pst(1, 30)= 87 - pst(2, 30)= 1 - pst(1, 31)= 88 - pst(2, 31)= 2 - pst(1, 32)= 90 - pst(2, 32)= 4 - pst(1, 33)= 94 - pst(2, 33)= 8 - pst(1, 34)= 0 - pst(2, 34)= 0 - pst(1, 35)= 0 - pst(2, 35)= 0 - pst(1, 36)= 0 - pst(2, 36)= 0 - pst(1, 37)= 0 - pst(2, 37)= 0 - pst(1, 38)= 0 - pst(2, 38)= 0 - pst(1, 39)= 0 - pst(2, 39)= 0 - pst(1, 40)= 0 - pst(2, 40)= 0 - pst(1, 41)= 0 - pst(2, 41)= 0 - pst(1, 42)= 0 - pst(2, 42)= 0 - pst(1, 43)= 0 - pst(2, 43)= 0 - pst(1, 44)= 0 - pst(2, 44)= 0 - pst(1, 45)= 0 - pst(2, 45)= 0 - pst(1, 46)= 0 - pst(2, 46)= 0 - pst(1, 47)= 0 - pst(2, 47)= 0 - pst(1, 48)= 0 - pst(2, 48)= 0 - pst(1, 49)= 0 - pst(2, 49)= 0 - pst(1, 50)= 0 - pst(2, 50)= 0 - pst(1, 51)= 0 - pst(2, 51)= 0 - pst(1, 52)= 0 - pst(2, 52)= 0 - pst(1, 53)= 0 - pst(2, 53)= 0 - pst(1, 54)= 0 - pst(2, 54)= 0 - pst(1, 55)= 0 - pst(2, 55)= 0 - pst(1, 56)= 0 - pst(2, 56)= 0 - pst(1, 57)= 0 - pst(2, 57)= 0 - pst(1, 58)= 0 - pst(2, 58)= 0 - pst(1, 59)= 0 - pst(2, 59)= 0 - pst(1, 60)= 0 - pst(2, 60)= 0 - pst(1, 61)= 0 - pst(2, 61)= 0 - pst(1, 62)= 0 - pst(2, 62)= 0 - pst(1, 63)= 0 - pst(2, 63)= 0 - pst(1, 64)= 0 - pst(2, 64)= 0 - pst(1, 65)= 0 - pst(2, 65)= 0 - pst(1, 66)= 0 - pst(2, 66)= 0 - pst(1, 67)= 0 - pst(2, 67)= 0 - pst(1, 68)= 0 - pst(2, 68)= 0 - pst(1, 69)= 0 - pst(2, 69)= 0 - pst(1, 70)= 0 - pst(2, 70)= 0 - pst(1, 71)= 0 - pst(2, 71)= 0 - pst(1, 72)= 0 - pst(2, 72)= 0 - pst(1, 73)= 0 - pst(2, 73)= 0 - pst(1, 74)= 0 - pst(2, 74)= 0 - pst(1, 75)= 0 - pst(2, 75)= 0 - pst(1, 76)= 0 - pst(2, 76)= 0 - pst(1, 77)= 0 - pst(2, 77)= 0 - pst(1, 78)= 0 - pst(2, 78)= 0 - pst(1, 79)= 0 - pst(2, 79)= 0 - pst(1, 80)= 0 - pst(2, 80)= 0 - pst(1, 81)= 0 - pst(2, 81)= 0 - pst(1, 82)= 0 - pst(2, 82)= 0 - pst(1, 83)= 0 - pst(2, 83)= 0 - pst(1, 84)= 0 - pst(2, 84)= 0 - pst(1, 85)= 0 - pst(2, 85)= 0 - pst(1, 86)= 0 - pst(2, 86)= 0 - pst(1, 87)= 0 - pst(2, 87)= 0 - pst(1, 88)= 0 - pst(2, 88)= 0 - pst(1, 89)= 0 - pst(2, 89)= 0 - pst(1, 90)= 0 - pst(2, 90)= 0 - pst(1, 91)= 0 - pst(2, 91)= 0 - pst(1, 92)= 0 - pst(2, 92)= 0 - pst(1, 93)= 0 - pst(2, 93)= 0 - pst(1, 94)= 0 - pst(2, 94)= 0 - pst(1, 95)= 0 - pst(2, 95)= 0 - pst(1, 96)= 0 - pst(2, 96)= 0 - pst(1, 97)= 0 - pst(2, 97)= 0 - pst(1, 98)= 0 - pst(2, 98)= 0 - pst(1, 99)= 0 - pst(2, 99)= 0 - pst(1,100)= 0 - pst(2,100)= 0 - pst(1,101)= 0 - pst(2,101)= 0 - pst(1,102)= 0 - pst(2,102)= 0 - pst(1,103)= 0 - pst(2,103)= 0 - pst(1,104)= 0 - pst(2,104)= 0 - pst(1,105)= 0 - pst(2,105)= 0 - pst(1,106)= 0 - pst(2,106)= 0 - pst(1,107)= 0 - pst(2,107)= 0 - pst(1,108)= 0 - pst(2,108)= 0 - pst(1,109)= 0 - pst(2,109)= 0 - pst(1,110)= 0 - pst(2,110)= 0 - pst(1,111)= 0 - pst(2,111)= 0 - pst(1,112)= 0 - pst(2,112)= 0 - pst(1,113)= 0 - pst(2,113)= 0 - pst(1,114)= 0 - pst(2,114)= 0 - pst(1,115)= 0 - pst(2,115)= 0 - pst(1,116)= 0 - pst(2,116)= 0 - pst(1,117)= 0 - pst(2,117)= 0 - pst(1,118)= 0 - pst(2,118)= 0 - pst(1,119)= 0 - pst(2,119)= 0 - pst(1,120)= 0 - pst(2,120)= 0 - pst(1,121)= 0 - pst(2,121)= 0 - pst(1,122)= 0 - pst(2,122)= 0 - pst(1,123)= 0 - pst(2,123)= 0 - pst(1,124)= 0 - pst(2,124)= 0 - pst(1,125)= 0 - pst(2,125)= 0 - pst(1,126)= 0 - pst(2,126)= 0 - pst(1,127)= 0 - pst(2,127)= 0 - pst(1,128)= 0 - pst(2,128)= 0 - pst(1,129)= 0 - pst(2,129)= 0 - pst(1,130)= 0 - pst(2,130)= 0 - pst(1,131)= 0 - pst(2,131)= 0 - pst(1,132)= 0 - pst(2,132)= 0 - pst(1,133)= 0 - pst(2,133)= 0 - pst(1,134)= 0 - pst(2,134)= 0 - pst(1,135)= 0 - pst(2,135)= 0 - pst(1,136)= 0 - pst(2,136)= 0 - pst(1,137)= 0 - pst(2,137)= 0 - pst(1,138)= 0 - pst(2,138)= 0 - pst(1,139)= 0 - pst(2,139)= 0 - pst(1,140)= 0 - pst(2,140)= 0 - pst(1,141)= 0 - pst(2,141)= 0 - pst(1,142)= 0 - pst(2,142)= 0 - pst(1,143)= 0 - pst(2,143)= 0 - pst(1,144)= 0 - pst(2,144)= 0 - pst(1,145)= 0 - pst(2,145)= 0 - pst(1,146)= 0 - pst(2,146)= 0 - pst(1,147)= 0 - pst(2,147)= 0 - pst(1,148)= 0 - pst(2,148)= 0 - pst(1,149)= 0 - pst(2,149)= 0 - pst(1,150)= 0 - pst(2,150)= 0 - pst(1,151)= 0 - pst(2,151)= 0 - pst(1,152)= 0 - pst(2,152)= 0 - pst(1,153)= 0 - pst(2,153)= 0 - pst(1,154)= 0 - pst(2,154)= 0 - pst(1,155)= 0 - pst(2,155)= 0 - pst(1,156)= 0 - pst(2,156)= 0 - pst(1,157)= 0 - pst(2,157)= 0 - pst(1,158)= 0 - pst(2,158)= 0 - pst(1,159)= 0 - pst(2,159)= 0 - pst(1,160)= 0 - pst(2,160)= 0 - pst(1,161)= 0 - pst(2,161)= 0 - pst(1,162)= 0 - pst(2,162)= 0 - pst(1,163)= 0 - pst(2,163)= 0 - pst(1,164)= 0 - pst(2,164)= 0 - pst(1,165)= 0 - pst(2,165)= 0 - pst(1,166)= 0 - pst(2,166)= 0 - pst(1,167)= 0 - pst(2,167)= 0 - pst(1,168)= 0 - pst(2,168)= 0 - pst(1,169)= 0 - pst(2,169)= 0 - pst(1,170)= 0 - pst(2,170)= 0 - pst(1,171)= 0 - pst(2,171)= 0 - pst(1,172)= 0 - pst(2,172)= 0 - pst(1,173)= 0 - pst(2,173)= 0 - pst(1,174)= 0 - pst(2,174)= 0 - pst(1,175)= 0 - pst(2,175)= 0 - pst(1,176)= 0 - pst(2,176)= 0 - pst(1,177)= 0 - pst(2,177)= 0 - pst(1,178)= 0 - pst(2,178)= 0 - pst(1,179)= 0 - pst(2,179)= 0 - pst(1,180)= 0 - pst(2,180)= 0 - pst(1,181)= 0 - pst(2,181)= 0 - pst(1,182)= 0 - pst(2,182)= 0 - pst(1,183)= 0 - pst(2,183)= 0 - pst(1,184)= 0 - pst(2,184)= 0 - pst(1,185)= 0 - pst(2,185)= 0 - pst(1,186)= 0 - pst(2,186)= 0 - pst(1,187)= 0 - pst(2,187)= 0 - pst(1,188)= 0 - pst(2,188)= 0 - pst(1,189)= 0 - pst(2,189)= 0 - pst(1,190)= 0 - pst(2,190)= 0 - pst(1,191)= 0 - pst(2,191)= 0 - pst(1,192)= 0 - pst(2,192)= 0 - pst(1,193)= 0 - pst(2,193)= 0 - pst(1,194)= 0 - pst(2,194)= 0 - pst(1,195)= 0 - pst(2,195)= 0 - pst(1,196)= 0 - pst(2,196)= 0 - pst(1,197)= 0 - pst(2,197)= 0 - pst(1,198)= 0 - pst(2,198)= 0 - pst(1,199)= 0 - pst(2,199)= 0 - pst(1,200)= 0 - pst(2,200)= 0 - pst(1,201)= 0 - pst(2,201)= 0 - pst(1,202)= 0 - pst(2,202)= 0 - pst(1,203)= 0 - pst(2,203)= 0 - pst(1,204)= 0 - pst(2,204)= 0 - pst(1,205)= 0 - pst(2,205)= 0 - pst(1,206)= 0 - pst(2,206)= 0 - pst(1,207)= 0 - pst(2,207)= 0 - pst(1,208)= 0 - pst(2,208)= 0 - pst(1,209)= 0 - pst(2,209)= 0 - pst(1,210)= 0 - pst(2,210)= 0 - pst(1,211)= 0 - pst(2,211)= 0 - pst(1,212)= 0 - pst(2,212)= 0 - pst(1,213)= 0 - pst(2,213)= 0 - pst(1,214)= 0 - pst(2,214)= 0 - pst(1,215)= 0 - pst(2,215)= 0 - pst(1,216)= 0 - pst(2,216)= 0 - pst(1,217)= 0 - pst(2,217)= 0 - pst(1,218)= 0 - pst(2,218)= 0 - pst(1,219)= 0 - pst(2,219)= 0 - pst(1,220)= 0 - pst(2,220)= 0 - pst(1,221)= 0 - pst(2,221)= 0 - pst(1,222)= 0 - pst(2,222)= 0 - pst(1,223)= 0 - pst(2,223)= 0 - pst(1,224)= 0 - pst(2,224)= 0 - pst(1,225)= 0 - pst(2,225)= 0 - pst(1,226)= 0 - pst(2,226)= 0 - pst(1,227)= 0 - pst(2,227)= 0 - pst(1,228)= 0 - pst(2,228)= 0 - pst(1,229)= 0 - pst(2,229)= 0 - pst(1,230)= 0 - pst(2,230)= 0 - pst(1,231)= 0 - pst(2,231)= 0 - pst(1,232)= 0 - pst(2,232)= 0 - pst(1,233)= 0 - pst(2,233)= 0 - pst(1,234)= 0 - pst(2,234)= 0 - pst(1,235)= 0 - pst(2,235)= 0 - pst(1,236)= 0 - pst(2,236)= 0 - pst(1,237)= 0 - pst(2,237)= 0 - pst(1,238)= 0 - pst(2,238)= 0 - pst(1,239)= 0 - pst(2,239)= 0 - pst(1,240)= 0 - pst(2,240)= 0 - pst(1,241)= 0 - pst(2,241)= 0 - pst(1,242)= 0 - pst(2,242)= 0 - pst(1,243)= 0 - pst(2,243)= 0 - pst(1,244)= 0 - pst(2,244)= 0 - pst(1,245)= 0 - pst(2,245)= 0 - pst(1,246)= 0 - pst(2,246)= 0 - pst(1,247)= 0 - pst(2,247)= 0 - pst(1,248)= 0 - pst(2,248)= 0 - pst(1,249)= 0 - pst(2,249)= 0 - pst(1,250)= 0 - pst(2,250)= 0 - pst(1,251)= 0 - pst(2,251)= 0 - pst(1,252)= 0 - pst(2,252)= 0 - pst(1,253)= 0 - pst(2,253)= 0 - pst(1,254)= 0 - pst(2,254)= 0 - pst(1,255)= 0 - pst(2,255)= 0 - pst(1,256)= 0 - pst(2,256)= 0 - pst(1,257)= 0 - pst(2,257)= 0 - pst(1,258)= 0 - pst(2,258)= 0 - pst(1,259)= 0 - pst(2,259)= 0 - pst(1,260)= 0 - pst(2,260)= 0 - pst(1,261)= 0 - pst(2,261)= 0 - pst(1,262)= 0 - pst(2,262)= 0 - pst(1,263)= 0 - pst(2,263)= 0 - pst(1,264)= 0 - pst(2,264)= 0 - pst(1,265)= 0 - pst(2,265)= 0 - pst(1,266)= 0 - pst(2,266)= 0 - pst(1,267)= 0 - pst(2,267)= 0 - pst(1,268)= 0 - pst(2,268)= 0 - pst(1,269)= 0 - pst(2,269)= 0 - pst(1,270)= 0 - pst(2,270)= 0 - pst(1,271)= 0 - pst(2,271)= 0 - pst(1,272)= 0 - pst(2,272)= 0 - pst(1,273)= 0 - pst(2,273)= 0 - pst(1,274)= 0 - pst(2,274)= 0 - pst(1,275)= 0 - pst(2,275)= 0 - pst(1,276)= 0 - pst(2,276)= 0 - pst(1,277)= 0 - pst(2,277)= 0 - pst(1,278)= 0 - pst(2,278)= 0 - pst(1,279)= 0 - pst(2,279)= 0 - pst(1,280)= 0 - pst(2,280)= 0 - pst(1,281)= 0 - pst(2,281)= 0 - pst(1,282)= 0 - pst(2,282)= 0 - pst(1,283)= 0 - pst(2,283)= 0 - pst(1,284)= 0 - pst(2,284)= 0 - pst(1,285)= 0 - pst(2,285)= 0 - pst(1,286)= 0 - pst(2,286)= 0 - pst(1,287)= 0 - pst(2,287)= 0 - pst(1,288)= 0 - pst(2,288)= 0 - pst(1,289)= 0 - pst(2,289)= 0 - pst(1,290)= 0 - pst(2,290)= 0 - pst(1,291)= 0 - pst(2,291)= 0 - pst(1,292)= 0 - pst(2,292)= 0 - pst(1,293)= 0 - pst(2,293)= 0 - pst(1,294)= 0 - pst(2,294)= 0 - pst(1,295)= 0 - pst(2,295)= 0 - pst(1,296)= 0 - pst(2,296)= 0 - pst(1,297)= 0 - pst(2,297)= 0 - pst(1,298)= 0 - pst(2,298)= 0 - pst(1,299)= 0 - pst(2,299)= 0 - pst(1,300)= 0 - pst(2,300)= 0 - pst(1,301)= 0 - pst(2,301)= 0 - pst(1,302)= 0 - pst(2,302)= 0 - pst(1,303)= 0 - pst(2,303)= 0 - pst(1,304)= 0 - pst(2,304)= 0 - pst(1,305)= 0 - pst(2,305)= 0 - pst(1,306)= 0 - pst(2,306)= 0 - pst(1,307)= 0 - pst(2,307)= 0 - pst(1,308)= 0 - pst(2,308)= 0 - pst(1,309)= 0 - pst(2,309)= 0 - pst(1,310)= 0 - pst(2,310)= 0 - pst(1,311)= 0 - pst(2,311)= 0 - pst(1,312)= 0 - pst(2,312)= 0 - pst(1,313)= 0 - pst(2,313)= 0 - pst(1,314)= 0 - pst(2,314)= 0 - pst(1,315)= 0 - pst(2,315)= 0 - pst(1,316)= 0 - pst(2,316)= 0 - pst(1,317)= 0 - pst(2,317)= 0 - pst(1,318)= 0 - pst(2,318)= 0 - pst(1,319)= 0 - pst(2,319)= 0 - pst(1,320)= 0 - pst(2,320)= 0 - pst(1,321)= 0 - pst(2,321)= 0 - pst(1,322)= 0 - pst(2,322)= 0 - pst(1,323)= 0 - pst(2,323)= 0 - pst(1,324)= 0 - pst(2,324)= 0 - pst(1,325)= 0 - pst(2,325)= 0 - pst(1,326)= 0 - pst(2,326)= 0 - pst(1,327)= 0 - pst(2,327)= 0 - pst(1,328)= 0 - pst(2,328)= 0 - pst(1,329)= 0 - pst(2,329)= 0 - pst(1,330)= 0 - pst(2,330)= 0 - pst(1,331)= 0 - pst(2,331)= 0 - pst(1,332)= 0 - pst(2,332)= 0 - pst(1,333)= 0 - pst(2,333)= 0 - pst(1,334)= 0 - pst(2,334)= 0 - pst(1,335)= 0 - pst(2,335)= 0 - pst(1,336)= 0 - pst(2,336)= 0 - pst(1,337)= 0 - pst(2,337)= 0 - pst(1,338)= 0 - pst(2,338)= 0 - pst(1,339)= 0 - pst(2,339)= 0 - pst(1,340)= 0 - pst(2,340)= 0 - pst(1,341)= 0 - pst(2,341)= 0 - pst(1,342)= 0 - pst(2,342)= 0 - pst(1,343)= 0 - pst(2,343)= 0 - pst(1,344)= 0 - pst(2,344)= 0 - pst(1,345)= 0 - pst(2,345)= 0 - pst(1,346)= 0 - pst(2,346)= 0 - pst(1,347)= 0 - pst(2,347)= 0 - pst(1,348)= 0 - pst(2,348)= 0 - pst(1,349)= 0 - pst(2,349)= 0 - pst(1,350)= 0 - pst(2,350)= 0 - pst(1,351)= 0 - pst(2,351)= 0 - pst(1,352)= 0 - pst(2,352)= 0 - pst(1,353)= 0 - pst(2,353)= 0 - pst(1,354)= 0 - pst(2,354)= 0 - pst(1,355)= 0 - pst(2,355)= 0 - pst(1,356)= 0 - pst(2,356)= 0 - pst(1,357)= 0 - pst(2,357)= 0 - pst(1,358)= 0 - pst(2,358)= 0 - pst(1,359)= 0 - pst(2,359)= 0 - pst(1,360)= 0 - pst(2,360)= 0 - pst(1,361)= 0 - pst(2,361)= 0 - pst(1,362)= 0 - pst(2,362)= 0 - pst(1,363)= 0 - pst(2,363)= 0 - pst(1,364)= 0 - pst(2,364)= 0 - pst(1,365)= 0 - pst(2,365)= 0 - pst(1,366)= 0 - pst(2,366)= 0 - pst(1,367)= 0 - pst(2,367)= 0 - pst(1,368)= 0 - pst(2,368)= 0 - pst(1,369)= 0 - pst(2,369)= 0 - pst(1,370)= 0 - pst(2,370)= 0 - pst(1,371)= 0 - pst(2,371)= 0 - pst(1,372)= 0 - pst(2,372)= 0 - pst(1,373)= 0 - pst(2,373)= 0 - pst(1,374)= 0 - pst(2,374)= 0 - pst(1,375)= 0 - pst(2,375)= 0 - pst(1,376)= 0 - pst(2,376)= 0 - pst(1,377)= 0 - pst(2,377)= 0 - pst(1,378)= 0 - pst(2,378)= 0 - pst(1,379)= 0 - pst(2,379)= 0 - pst(1,380)= 0 - pst(2,380)= 0 - pst(1,381)= 0 - pst(2,381)= 0 - pst(1,382)= 0 - pst(2,382)= 0 - pst(1,383)= 0 - pst(2,383)= 0 - pst(1,384)= 0 - pst(2,384)= 0 - pst(1,385)= 0 - pst(2,385)= 0 - pst(1,386)= 0 - pst(2,386)= 0 - pst(1,387)= 0 - pst(2,387)= 0 - pst(1,388)= 0 - pst(2,388)= 0 - pst(1,389)= 0 - pst(2,389)= 0 - pst(1,390)= 0 - pst(2,390)= 0 - pst(1,391)= 0 - pst(2,391)= 0 - pst(1,392)= 0 - pst(2,392)= 0 - pst(1,393)= 0 - pst(2,393)= 0 - pst(1,394)= 0 - pst(2,394)= 0 - pst(1,395)= 0 - pst(2,395)= 0 - pst(1,396)= 0 - pst(2,396)= 0 - pst(1,397)= 0 - pst(2,397)= 0 - pst(1,398)= 0 - pst(2,398)= 0 - pst(1,399)= 0 - pst(2,399)= 0 - pst(1,400)= 0 - pst(2,400)= 0 -End SUBROUTINE init_dip_planar_data -End Module dip_param diff --git a/src/model/invariants_nh3.f90 b/src/model/invariants_nh3.f90 deleted file mode 100644 index b42dddd..0000000 --- a/src/model/invariants_nh3.f90 +++ /dev/null @@ -1,172 +0,0 @@ -Module invariants_mod - implicit none - contains - !---------------------------------------------------- - subroutine invariants(a,xs,ys,xb,yb,inv) - implicit none - !include "nnparams.incl" - double precision, intent(in) :: a, xs, ys, xb, yb - double precision, intent(out) ::inv(3) - double precision:: invar(23) - !double precision:: invar(24) - complex(8) :: q1, q2 - LOGICAL,PARAMETER:: debg =.FALSE. - integer :: i - ! express the coordinate in complex - - q1 = dcmplx(xs, ys) - q2 = dcmplx(xb, yb) - - - ! compute the invariants - !inv( 4) = a - !inv( 4) = b**2 - - ! INVARIANTS OF KIND II - !------------------------ - - invar(1) = dreal( q1 * conjg(q1) ) ! r11 - invar(2) = dreal( q1 * conjg(q2) ) ! r12 - invar(3) = dreal( q2 * conjg(q2) ) ! r22 - invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2 - - - !INVATIANTS OF KIND III - !------------------------ - - invar(5) = dreal( q1 * q1 * q1 ) ! r111 - invar(6) = dreal( q1 * q1 * q2 ) ! r112 - invar(7) = dreal( q1 * q2 * q2 ) ! r122 - invar(8) = dreal( q2 * q2 * q2 ) ! r222 - invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2 - invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2 - invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2 - invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222 - - ! INVARIANTS OF KIND IV - !------------------------- - - invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 )) - invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 )) - invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 )) - invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 )) - - ! INVARIANTS OF KIND V - !---------------------- - - invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 )) - invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 )) - invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 )) - invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 )) - invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 )) - invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 )) - invar(23) = a - - ! the only non zero invariant for bend pure cuts - - inv(1) = invar( 1) - inv(2) = invar( 5) - inv(3) = invar( 9) - ! inv(4) = invar(1) - ! inv(5) = invar(5) - ! inv(6) = invar(9) - !inv(17)=invar(18) - !inv(19) = invar() - !inv(2) = invar(5) - !inv(3) = invar(9) - !inv(1:22)=invar(1:22) - !inv(4:8) = invar(5:9) - !inv(9)=invar(12) - - - if (debg) then - write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) - write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) - write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) - write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) - - write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES" - write(14,*)"---------------------------------------" - write(14,*)"xs =",dreal(q1), "ys=",dimag(q1) - endif - ! modify the invariants to only consider few of them - ! - !invar(13:22)=0.0d0 - end subroutine invariants - - -! subroutine invariants_ch(a, xs, ys, xb, yb, invar) -! implicit none -! double precision, intent(in) :: a, xs, ys, xb, yb -! double precision, intent(out) :: invar(25) -! -! complex(8) :: q1, q2 -! double precision :: r11, r22, r12, eta -! double precision :: I1, I2, I3, I4 -! -! ! Define complex coordinates -! q1 = dcmplx(xs, ys) -! q2 = dcmplx(xb, yb) -! -! ! Quadratic invariants -! r11 = dreal(q1*conjg(q1)) -! r22 = dreal(q2*conjg(q2)) -! r12 = dreal(q1*conjg(q2)) -! eta = dimag(q1*conjg(q2)) -! -! ! Cubic imag parts (used many times) -! I1 = dimag(q1*q1*q1) -! I2 = dimag(q1*q1*q2) -! I3 = dimag(q1*q2*q2) -! I4 = dimag(q2*q2*q2) -! ! the totally symmetric invariant 25 -! invar(25) = a -! ! Store invariants -! invar(1) = r11 -! invar(2) = r12 -! invar(3) = r22 -! invar(4) = eta*eta -! -! invar(5) = dreal(q1*q1*q1) -! invar(6) = dreal(q1*q1*q2) -! invar(7) = dreal(q1*q2*q2) -! invar(8) = dreal(q2*q2*q2) -! -! invar(9) = I1*I1 -! invar(10) = I2*I2 -! invar(11) = I3*I3 -! invar(12) = I4*I4 -! -! invar(13) = eta * I1 -! invar(14) = eta * I2 -! invar(15) = eta * I3 -! invar(16) = eta * I4 -! -! ! Pure quartic invariants -! invar(17) = r11*r11 + r22*r22 -! invar(18) = r11*r22 - r12*r12 - eta*eta -! -! ! Cubic-imag bilinear invariants (minimal independent set) -! invar(19) = I1*I2 -! invar(20) = I1*I3 -! invar(21) = I1*I4 -! invar(22) = I2*I3 -! invar(23) = I2*I4 -! -! ! Odd cubic invariant -! invar(24) = eta**3 -! -! if (debg) then -! write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) -! write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) -! write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) -! write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) -! -! write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES" -! write(14,*)"---------------------------------------" -! write(14,*)"xs =",dreal(q1), "ys=",dimag(q1) -! endif -! -! end subroutine - -end module invariants_mod diff --git a/src/model/newmodes.f b/src/model/newmodes.f deleted file mode 100644 index 1dfb0f7..0000000 --- a/src/model/newmodes.f +++ /dev/null @@ -1,87 +0,0 @@ -! Author: jnshuti -! Created: 2025-10-24 16:45:07 -! Last modified: 2025-10-24 16:45:27 jnshuti -! Modes: -! Mode 1 - mode( 1, 1) = 0.0000000000d0 - mode( 2, 1) = 0.0000000000d0 - mode( 3, 1) = 0.0000000000d0 - mode( 4, 1) = 0.5773502692d0 - mode( 5, 1) = 0.0000000000d0 - mode( 6, 1) = 0.0000000000d0 - mode( 7, 1) = -0.2886751346d0 - mode( 8, 1) = 0.5000000000d0 - mode( 9, 1) = 0.0000000000d0 - mode(10, 1) = -0.2886751346d0 - mode(11, 1) = -0.5000000000d0 - mode(12, 1) = 0.0000000000d0 - -! Mode 2 - mode( 1, 2) = 0.3121511560d0 - mode( 2, 2) = 0.0000000000d0 - mode( 3, 2) = 0.0000000000d0 - mode( 4, 2) = -0.7756982471d0 - mode( 5, 2) = 0.0000000000d0 - mode( 6, 2) = 0.0000000000d0 - mode( 7, 2) = -0.1939245618d0 - mode( 8, 2) = 0.3358871938d0 - mode( 9, 2) = 0.0000000000d0 - mode(10, 2) = -0.1939245618d0 - mode(11, 2) = -0.3358871938d0 - mode(12, 2) = 0.0000000000d0 - -! Mode 3 - mode( 1, 3) = 0.0000000000d0 - mode( 2, 3) = 0.3121511560d0 - mode( 3, 3) = 0.0000000000d0 - mode( 4, 3) = 0.0000000000d0 - mode( 5, 3) = 0.0000000000d0 - mode( 6, 3) = 0.0000000000d0 - mode( 7, 3) = 0.3358871938d0 - mode( 8, 3) = -0.5817736853d0 - mode( 9, 3) = 0.0000000000d0 - mode(10, 3) = -0.3358871938d0 - mode(11, 3) = -0.5817736853d0 - mode(12, 3) = 0.0000000000d0 - -! Mode 4 - mode( 1, 4) = 0.2830826954d0 - mode( 2, 4) = 0.0000000000d0 - mode( 3, 4) = 0.0000000000d0 - mode( 4, 4) = 0.0759441281d0 - mode( 5, 4) = 0.0000000000d0 - mode( 6, 4) = 0.0000000000d0 - mode( 7, 4) = -0.5655692225d0 - mode( 8, 4) = -0.3703779057d0 - mode( 9, 4) = 0.0000000000d0 - mode(10, 4) = -0.5655692225d0 - mode(11, 4) = 0.3703779057d0 - mode(12, 4) = 0.0000000000d0 - -! Mode 5 - mode( 1, 5) = 0.0000000000d0 - mode( 2, 5) = 0.2830826954d0 - mode( 3, 5) = 0.0000000000d0 - mode( 4, 5) = 0.0000000000d0 - mode( 5, 5) = -0.7794070061d0 - mode( 6, 5) = 0.0000000000d0 - mode( 7, 5) = -0.3703779057d0 - mode( 8, 5) = -0.1378936554d0 - mode( 9, 5) = 0.0000000000d0 - mode(10, 5) = 0.3703779057d0 - mode(11, 5) = -0.1378936554d0 - mode(12, 5) = 0.0000000000d0 - -! Mode 6 - mode( 1, 6) = 0.0000000000d0 - mode( 2, 6) = 0.0000000000d0 - mode( 3, 6) = 0.4213954872d0 - mode( 4, 6) = 0.0000000000d0 - mode( 5, 6) = 0.0000000000d0 - mode( 6, 6) = -0.5235856642d0 - mode( 7, 6) = 0.0000000000d0 - mode( 8, 6) = 0.0000000000d0 - mode( 9, 6) = -0.5235856642d0 - mode(10, 6) = 0.0000000000d0 - mode(11, 6) = 0.0000000000d0 - mode(12, 6) = -0.5235856642d0 diff --git a/src/model/nnadia_nh3.f b/src/model/nnadia_nh3.f deleted file mode 100644 index d9cc935..0000000 --- a/src/model/nnadia_nh3.f +++ /dev/null @@ -1,87 +0,0 @@ - subroutine nninit_mod() - use nn_params - implicit none - - end subroutine - - - subroutine nnadia(coords,coeffs,adiaoutp) - USE pesmod, only: potential_nh3_low - Use diabmodel, only: diab_x - - use nn_params - use nncommons - implicit none - ! returns THE DIABATIC MATRIX UPPER MATRIX OF BOTH X AND Y COMPONENT OF DIPOLE -! -! coords: vector containing symmetry-adapted coordinates. -! coeffs: vector containing coeffs of diab dipole matrix - -! adiaoutp: upper diab matrix of dipole -! x component come first and then y component -! -! dmat_x & dmat_y: analytic model of dipole, x and y -! -! matdim: dimension of hamiltoniam matrix -! mkeigen: integer: no eigenvectors if =0 -! eigenvs: dummy storage for eigenvectors - - !include 'JTmod.incl' - DOUBLE PRECISION,intent(in):: coords(maxnin) - DOUBLE PRECISION,INTENT(INOUT)::coeffs(maxnout) - DOUBLE PRECISION,INTENT(OUT):: adiaoutp(maxpout) - DOUBLE PRECISION, DIMENSION(ndiab,ndiab):: v,U,ex,temp - double precision, dimension(ndiab):: e - Double precision, allocatable:: mat(:,:) - integer i, ij, j, max_row - - ! variable for lapack - integer,parameter :: lwork = 1000 - double precision work(lwork) - integer info - - call potential_nh3_low(V,coords,coeffs) - - ! diagonalize the pes - u=0.0d0 - allocate(mat,source=V) - call DSYEV('V','U',ndiab,mat,ndiab,e,work,lwork, info) - - if (info .ne. 0) then - write(6,*)"The diagonalizing routine failed", info - stop - end if - U(:,:)= mat(:,:) - deallocate(mat) - do i =1,ndiab - adiaoutp(i)=e(i) - enddo - - ! Fix the phase of U matrix - ! fix phas of U - do i =1,ndiab - max_row=maxloc(abs(U(:,i)),1) - if(U(max_row,i) .lt. 0) then - U(:,i) = -1.0d0*U(:,i) - endif - enddo - - call diab_x(ex,coords,coeffs) ! ex is the dipole x component - - ! transforrm the dipole - - temp = matmul(transpose(U),ex) - ex = matmul(temp,U) - ij = 4 - do i =1,ndiab - do j = i,ndiab - ij = ij+1 - adiaoutp(ij) = ex(i,j) - !adiaoutp(ij) = 0.0d0 - enddo - enddo - - END SUBROUTINE - - - diff --git a/src/model/nncoords_nh3.f90 b/src/model/nncoords_nh3.f90 deleted file mode 100644 index f70883e..0000000 --- a/src/model/nncoords_nh3.f90 +++ /dev/null @@ -1,91 +0,0 @@ - module trans_coord - implicit none - contains - - subroutine cart2mode(qq) - use invariants_mod, only: invariants - use nn_params - use nncommons - implicit none - double precision, intent(inout):: qq(maxnin) - double precision :: q(12),inv(3),qm(6) - double precision, parameter :: m_N=25526.04237518618d0 - double precision, parameter :: m_H=1837.152647439619d0 - double precision, parameter :: Ang2Bohr=1.889726124565d0 - - double precision diff(12), ref_geom(12) - double precision q_mod(12) -! transformation matrix - double precision mode(12,6) - integer i -! final internal coord vector - - - include 'newmodes.f' - - ref_geom(1:12)=0.d0 - ref_geom( 4 )=1.022871078d0 - ref_geom( 7 )=-.5d0*ref_geom(4) - !ref_geom( 8 )=sqrt(1.5d0)*ref_geom(4) - ref_geom( 8)=0.5d0*sqrt(3.0d0)*ref_geom(4) - ref_geom(10 )=-.5d0*ref_geom(4) - !ref_geom(11 )=-sqrt(1.5d0)*ref_geom(4) - ref_geom(11)=-0.5d0*sqrt(3.0d0)*ref_geom(4) - - !massN=25526.04237518618 - !massH=1837.152647439619 - ! - q(1:12)=qq(1:12) - q_mod = (1.0d0/Ang2Bohr)*q - -! difference vector - diff=q_mod-ref_geom - -! mass weighting - diff(1:3)=diff(1:3)*sqrt(m_N) - diff(4:12)=diff(4:12)*sqrt(m_H) - - qm=matmul(transpose(mode),diff) - do i =1,size(qm) - if (abs(qm(i)) .lt. 1.0d-6) qm(i)=0.0d0 - enddo - !qm(3) = -qm(3) - call invariants(qm(1),qm(2),qm(3),qm(4),qm(5),inv) - qq =0.0d0 - qq(1:len_in)=inv(1:len_in) - qq(len_in+1:len_in+6)=qm(1:6) - - - end subroutine cart2mode - end module trans_coord -!**** Define coordinate transformation applied to the input before fit. -!*** -!*** -!****Conventions: -!*** -!*** ctrans: subroutine transforming a single point in coordinate space - - subroutine trans_in(pat_in) - !use ctrans_mod, only: ctrans - use trans_coord, only: cart2mode - use nn_params - use nncommons - implicit none - - !include 'nndbg.incl' - - double precision pat_in(maxnin,maxpats) - !integer ntot - - integer j,i,jj - jj=1 - do j=1,sets - !write(62,*)"Scan ",j - do i =1,ndata(j) - call cart2mode(pat_in(:,jj)) - ! FIRST ELEMENT OF PAT-IN ARE USED BY NEURON NETWORK - !write(62,'(*(ES16.8))') pat_in(:,jj) - jj = jj +1 - enddo - enddo - end diff --git a/src/model/pes_model.f90 b/src/model/pes_model.f90 deleted file mode 100644 index be53d68..0000000 --- a/src/model/pes_model.f90 +++ /dev/null @@ -1,172 +0,0 @@ -module pesmod - !use dim_parameter,only:qn,ndiab,pst - use iso_fortran_env, only: idp => int32, dp => real64 - use dip_param - - use nn_params - use nncommons, only: len_in - implicit none - - integer(idp),parameter:: ndiab=4 - private - public :: potential_nh3_low - contains - - - subroutine potential_nh3_low(e,q,nn_out) - implicit none - real(dp),intent(in)::q(maxnin) - real(dp),intent(out)::e(:,:) - real(dp),intent(inout):: nn_out(maxnout) - !logical,intent(in):: sym - !integer(idp),intent(in)::key - integer(idp) id,i, ii - real(dp) tmp,xs,xb,ys,yb,a,b,ss,sb - logical, parameter:: switch = .false. ! parameter to switch when you want to know Ex and Ey - logical, parameter:: sym = .false. ! to include r**2 in the Genetic fitting - real(dp),parameter:: tol = 1.0d-9 - real(dp),dimension(16):: scal, shift - !real(dp),dimension(17):: shi,scalee - !include "shift-scale-70N.incl" - !include "sh-scal-50.incl" - xs=q(len_in+2) - ys=q(len_in+3) - xb=q(len_in+4) - yb=q(len_in+5) - a=q(len_in+1) - b=q(len_in+6) - - ss= xs**2+ys**2 - sb= xb**2+yb**2 - call init_dip_planar_data() - shift =1.0d0 - scal = 1.0d-3 - !scal = 1.0d-3 - !scal(1:3) = 1.0d-3 - ii = 1 - do i =1 ,27 ! pes parameter key index - if( abs(p(i)) .gt. tol ) then - p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii)) - ii = ii + 1 - - endif - enddo - !write(6,*)" PES non zero", ii-1 - - - e = 0.0_dp - id = 1 ! 1 - ! v -term - ! oth order only - ! A2" - e(1,1) = e(1,1) + p(pst(1,id)) !- p(pst(1,1)) - if (sym) then - e(1,1)=e(1,1) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*sb - endif - id = id +1 !2 - ! E block - e(2,2) = e(2,2)+ p(pst(1,id))! - p(pst(1,1)) - e(3,3) = e(3,3)+ p(pst(1,id))! - p(pst(1,1)) - if (sym) then - e(2,2)=e(2,2)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 - e(3,3)=e(3,3)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 - endif - - id = id +1 ! 3 - - !A1' - e(4,4) = e(4,4) + p(pst(1,id))! - p(pst(1,1)) - if (sym) then - e(4,4)=e(4,4) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 & - + p(pst(1,id)+3)*ss**3 - endif - - ! W and z of JT block - - ! first order - id = id +1 !4 - - 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 ! 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) - e(2,3) = e(2,3) - p(pst(1,id))*(2*xs*ys) & - - p(pst(1,id)+1)*(2*xb*yb) & - - p(pst(1,id)+2)*(xs*yb+xb*ys) - - ! the coupling of A2' and E' - ! order 1 - id = id +1 !6 - - e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*yb) - e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb) - ! order 2 - id =id +1 ! 7 - e(1,2) = e(1,2) +b*(p(pst(1,id))*(xs**2-ys**2) & - + p(pst(1,id)+1)*(xb**2-yb**2) & - + p(pst(1,id)+2)*(xs*xb -ys*yb)) - e(1,3) = e(1,3) +b*(p(pst(1,id))*(2*xs*ys) & - + p(pst(1,id)+1)*(2*xb*yb) & - + p(pst(1,id)+2)*(xs*yb + xb*ys)) - - ! the coupling of A1' and E' - id = id +1 ! 8 - e(2,4) = e(2,4) + p(pst(1,id))*xs + p(pst(1,id)+1)*xb - e(3,4) = e(3,4) - p(pst(1,id))*ys - p(pst(1,id)+1)*yb - if (sym) then - e(2,4) = e(2,4)+ p(pst(1,id)+2)*xs*ss - e(3,4) = e(3,4)- p(pst(1,id)+2)*ys*ss - endif - - ! order 2 - id = id +1 ! 9 - e(2,4) = e(2,4) + p(pst(1,id))*(xs**2-ys**2) & - + p(pst(1,id)+1)*(xb**2-yb**2) & - + p(pst(1,id)+2)*(xs*xb -ys*yb) - e(3,4) = e(3,4) + p(pst(1,id))*(2*xs*ys) & - + p(pst(1,id)+1)*(2*xb*yb) & - + p(pst(1,id)+2)*(xs*yb + xb*ys) - ! the couplign between A - id = id +1 ! 10 - e(1,4) = e(1,4) + b*p(pst(1,id)) - !write(*,*) 'Id = ',id - - call copy_2_lower_triangle(e) - if (switch) then - write( 6,*)" EX AND EY ARE SWITCHED" - tmp = e(2,2) - e(2,2) = e(3,3) - e(3,3) = tmp - - ! the coupling - - tmp = e(2,4) - e(2,4) = e(3,4) - e(3,4) = tmp - e(4,2) = e(2,4) - e(4,3) = e(3,4) - endif - end subroutine potential_nh3_low - - - subroutine copy_2_lower_triangle(mat) - real(dp), intent(inout) :: mat(:, :) - integer :: m, n - ! write lower triangle of matrix symmetrical - do n = 2, size(mat, 1) - do m = 1, n - 1 - mat(n, m) = mat(m, n) - end do - end do - end subroutine copy_2_lower_triangle - - -end module diff --git a/src/model/red_invariants.f90 b/src/model/red_invariants.f90 deleted file mode 100644 index 278051f..0000000 --- a/src/model/red_invariants.f90 +++ /dev/null @@ -1,72 +0,0 @@ -Module invariants_mod - implicit none - contains - !---------------------------------------------------- - subroutine invariants(a,xs,ys,xb,yb,b,invar) - implicit none - double precision, intent(in) :: a, xs, ys, xb, yb, b - double precision, intent(out) :: invar(24) - complex(8) :: q1, q2 - LOGICAL,PARAMETER:: debg =.FALSE. - integer :: i - ! express the coordinate in complex - - q1 = dcmplx(xs, ys) - q2 = dcmplx(xb, yb) - - ! compute the invariants - invar(24) = a - invar(23) =b**2 - - ! INVARIANTS OF KIND II - !------------------------ - - invar(1) = dreal( q1 * conjg(q1) ) ! r11 - invar(2) = dreal( q1 * conjg(q2) ) ! r12 - invar(3) = dreal( q2 * conjg(q2) ) ! r22 - invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2 - - - !INVATIANTS OF KIND III - !------------------------ - - invar(5) = dreal( q1 * q1 * q1 ) ! r111 - invar(6) = dreal( q1 * q1 * q2 ) ! r112 - invar(7) = dreal( q1 * q2 * q2 ) ! r122 - invar(8) = dreal( q2 * q2 * q2 ) ! r222 - invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2 - invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2 - invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2 - invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222 - - ! INVARIANTS OF KIND IV - !------------------------- - - invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 )) - invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 )) - invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 )) - invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 )) - - ! INVARIANTS OF KIND V - !---------------------- - - invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 )) - invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 )) - invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 )) - invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 )) - invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 )) - invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 )) - - if (debg) then - write(*,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4) - write(*,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12) - write(*,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16) - write(*,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22) - - write(*,*)"THE INPUT COORDINATE IN COMPLEX REPRES" - write(*,*)"---------------------------------------" - write(*,*)"xs =",dreal(q1), "ys=",dimag(q1) - endif - - end subroutine invariants -end module invariants_mod diff --git a/src/model/sh-scal-50.incl b/src/model/sh-scal-50.incl deleted file mode 100644 index 8001743..0000000 --- a/src/model/sh-scal-50.incl +++ /dev/null @@ -1,9 +0,0 @@ - -! shift and scale for 70N structrure -! Real(dp),parameter:: scal(17) = 1.0d0 ,shift(17) = 1.0d0 - Real(dp),parameter:: scal(17) = & - [0.00434,0.00239,0.00201,4.13538,9.04704,1.07732,0.79253, & - 1.86682,1.87266,0.58857,0.90899,2.07755,1.81564,1.31957,1.68503,1.48391,1.79901] - Real(dp),parameter:: shift(17) = & - [-0.00016,-0.00020, 0.00008,-0.04261,13.00281,-0.44002,-0.18611, & - 0.73999, 0.34618,-0.44183,-0.06997, 0.39988,-1.42376, 0.77752,-1.23819,-0.22926, 0.82597] diff --git a/src/model/shift-scale-70N.incl b/src/model/shift-scale-70N.incl deleted file mode 100644 index 31d948e..0000000 --- a/src/model/shift-scale-70N.incl +++ /dev/null @@ -1,11 +0,0 @@ - -! shift and scale for 70N structrure -! Real(dp),parameter:: scal(17) = 1.0d0 ,shift(17) = 1.0d0 - Real(dp),parameter:: scal(17) = & - [0.00401, 0.00265, 0.00335, 1.19496, 3.79876, 0.32139, 0.46267, & - 1.36249, 0.90118, 0.94406, 0.29014, 0.72172, 1.25435, 2.09880, & - 0.93696, 1.69131, 0.75050 ] - Real(dp),parameter:: shift(17) = & - [-0.00018,-0.00018,-0.00002,-0.03463,-0.74584, 0.09806, 0.16350, & - -0.96733,-0.62850, 0.57182,-0.00110, 0.10584, 0.29949,-0.82418, & - 0.57450, 0.04471, 0.28142 ] diff --git a/src/model/write_error.f90 b/src/model/write_error.f90 deleted file mode 100644 index d8f2ddd..0000000 --- a/src/model/write_error.f90 +++ /dev/null @@ -1,235 +0,0 @@ -module print_error - use iso_fortran_env, only: idp => int32, dp => real64 - use nn_params - use nncommons - implicit none - - -contains - - !======================================================================= - ! Compute and print error summary (RMS) for PES energies and dipole moments - ! Separately for total, weighted, energies (1:4), dipoles (5:inp_out) - !======================================================================= - subroutine print_ErrorSummary(par, pat_in, pat_out,ref_in,ref_out, & - wterr, typop, laystr, weistr, nlay, & - npat, nref) - - real(dp), intent(in) :: pat_in(maxnin, maxpats) - real(dp), intent(in) :: pat_out(maxpout, maxpats) - real(dp), intent(in) :: par(wbcap, maxset) - real(dp), intent(in) :: wterr(maxpout, maxpats) - real(dp), intent(in) :: ref_in(maxnin,maxpats), ref_out(maxpout,maxpats) - integer(idp), intent(in) :: laystr(3, maxlay) - integer(idp), intent(in) :: typop(maxtypes, maxlay) - integer(idp), intent(in) :: weistr(2, maxlay, 2) - integer(idp), intent(in) :: nlay, npat, nref - - ! Local variables - integer(idp) :: pt, k,i - - integer(idp) :: ntot - integer(idp),dimension(maxpout) :: cnt_tot,cnt_fit,cnt_val - real(dp), dimension(maxpout) :: total_rms_fit, total_rms_val,wrms_fit,wrms_val - real(dp) :: nnoutp(maxnout) - real(dp) :: ymod(maxpout, maxpats) - real(dp) :: wsum_fit - real(dp) :: rms_total, rms_weight - real(dp) :: en_rms, dip_rms - Real(dp) :: rms_total_fit, rms_w_fit - Real(dp) :: rms_total_val, rms_w_val - Real(dp) :: en_rms_fit,en_rms_val,dip_rms_fit,dip_rms_val - integer(idp), parameter :: std_out = 6 - - ! Total number of points (fitting + validation/reference) - - - ntot = npat + nref - - ! keep in mind that npat is the total number of point fit + ref - - !write(6,*)"Npat in the error summary", npat,nref - write(11,'(*(ES16.7))')par(:,1) - - ! Initialize accumulators - cnt_fit = 0 - cnt_tot = 0 - cnt_val = 0 - wsum_fit =0.0_dp - total_rms_fit(:) = 0.0_dp - wrms_fit(:) = 0.0_dp - total_rms_val (:) = 0.0_dp - wrms_val(:) = 0.0_dp - !write(11,'(*(ES18.6))')ref_in(1:3,:) - - ! Generate model predictions and accumulate squared errors - pt = 0 - - do pt =1, npat - if (pt > maxpats) then - write(std_out, *) "ERROR: pt exceeds maxpats!" - stop - end if - - ! Get network output - call neunet(pat_in(:, pt), nnoutp, par,typop, laystr, weistr, nlay) - - ! Transform to adiabatic representation (energies + dipoles) - call nnadia(pat_in(:, pt), nnoutp, ymod(:, pt)) - - ! Accumulate errors only for components with positive weight - do k = 1, inp_out - if (wterr(k, pt) > 0.0_dp) then - !cnt_tot(k) = cnt_tot(k) + 1 - ! Fitting / training set - cnt_fit(k) = cnt_fit(k) + 1 - total_rms_fit(k) = total_rms_fit(k) + (ymod(k,pt) - pat_out(k,pt))**2 - wrms_fit(k) = wrms_fit(k) + ((ymod(k,pt) - pat_out(k,pt))**2 & - * (wterr(k,pt))) - wsum_fit = wsum_fit + wterr(k,pt) - - - end if - end do - end do - ! validation set - do i =1, nref - !if (pt > maxpats) stop "Error pt exceeds maxpats" - pt = npat + i - call neunet(ref_in(:,i),nnoutp,par,typop,laystr,weistr,nlay) - - ! call the model - call nnadia(ref_in(:,i),nnoutp,ymod(:,pt)) - - ! calculate the errors - do k=1,inp_out - if (wterr(k,pt) > 0.0_dp) then - cnt_val (k) = cnt_val(k) + 1 - total_rms_val(k) = total_rms_val(k) + (ymod(k,pt) - ref_out(k,i))**2 - wrms_val(k) = wrms_val(k) + ((ymod(k,pt) - ref_out(k,i))**2 * & - (wterr(k,pt))) - !write(6,*)" IN here wter" - endif - enddo - !write(11,'(*(ES18.6))')wterr(:,pt) - - enddo - - - - ! Safety check: did we process the expected number of points? - - - ! Compute RMS values - ! TOTAL FIT +REF - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - rms_total = 0.0_dp - if ((sum(cnt_fit(1:inp_out)) + sum(cnt_val(1:inp_out))) > 0) then - rms_total = dsqrt( sum(total_rms_fit + total_rms_val) / & - real((sum(cnt_fit(1:inp_out))+sum(cnt_val(1:inp_out))), dp) ) - end if - - rms_weight = 0.0_dp - if (sum(wterr(1:inp_out, 1:npat+nref)) > 1.0e-12_dp) then ! avoid div-by-zero - rms_weight = dsqrt( sum(wrms_fit) + sum(wrms_val)) / sum(wterr(1:inp_out, 1:npat+nref)) - - !rms_weight = dsqrt(wt_rms /weighted_wt ) - end if - - en_rms = 0.0_dp - if ((sum(cnt_fit(1:4)) + sum(cnt_val(1:4))) > 0) then - - en_rms = sqrt( ( sum(total_rms_fit(1:4)) + sum(total_rms_val(1:4)) ) / & - real(sum(cnt_fit(1:4)) + sum(cnt_val(1:4)), dp) ) - - end if - - dip_rms = 0.0_dp - if ((sum(cnt_fit(5:inp_out)) + sum(cnt_val(5:inp_out)))> 0) then - dip_rms = sqrt( (sum(total_rms_fit(5:inp_out)) + sum(total_rms_val(5:inp_out))) / & - real((sum(cnt_fit(5:inp_out))+sum(cnt_val(5:inp_out))), dp) ) - end if - - rms_total_fit = 0.0_dp - if (sum(cnt_fit(1:inp_out)) > 0) then - rms_total_fit = dsqrt( sum(total_rms_fit(1:inp_out)) / real(sum(cnt_fit(1:inp_out)), dp) ) - endif - - rms_w_fit = 0.0_dp - if (wsum_fit > 1.0e-12_dp) then - rms_w_fit = sqrt( sum(wrms_fit(1:inp_out)) / wsum_fit ) - endif - - ! validation - rms_total_val = 0.0_dp - if (sum(cnt_val(1:inp_out)) > 0) then - rms_total_val = sqrt( sum(total_rms_val(1:inp_out)) / real(sum(cnt_val(1:inp_out)), dp) ) - endif - - rms_w_val = 0.0_dp - if (sum(wterr(1:inp_out,npat+1:npat+nref)) > 1.0e-12_dp) then - rms_w_val = sqrt( sum(wrms_val(1:inp_out)) / sum(wterr(1:inp_out,npat+1:npat+nref)) ) - endif - - ! compute each quantity error - en_rms_fit = merge( sqrt(sum(total_rms_fit(1:4))/real(sum(cnt_fit(1:4)),dp)),& - 0.0_dp, sum(cnt_fit(1:4))>0 ) - en_rms_val = merge( sqrt(sum(total_rms_val(1:4))/real(sum(cnt_val(1:4)),dp)), & - 0.0_dp, sum(cnt_val(1:4))>0 ) - - dip_rms_fit = merge( sqrt(sum(total_rms_fit(5:inp_out))/real(sum(cnt_fit(5:inp_out)),dp)), & - 0.0_dp, sum(cnt_fit(5:inp_out))>0 ) - dip_rms_val = merge( sqrt(sum(total_rms_val(5:inp_out))/real(sum(cnt_val(5:inp_out)),dp)), & - 0.0_dp, sum(cnt_val(5:inp_out))>0 ) - - ! Output - write(std_out, '(A)') " " - write(std_out, '(A)') "USER DEFINED ERROR METRIC" - write(std_out, '(A)') "######################################" - !write(std_out, '(A,ES18.8,A)') & - ! "Total RMS (all fitted components, unweighted) : ", rms_total * unit_con, & - ! " [ha]" - !write(std_out, '(A,ES18.8,A)') & - ! "Total Weighted RMS (fit + Ref) : ", rms_weight * unit_con - !write(std_out, '(A,ES18.8,A)') & - ! "Energy RMS : ", en_rms * unit_con, " [ha] " - !write(std_out, '(A,ES18.8,A)') & - ! "Dipole RMS : ", dip_rms * unit_con, " [ha] " - - ! compute the metric for fit data - - write(std_out, '(A)') "FITTING / TRAINING SET" - write(std_out, '(A)') "#############################" - write(std_out,'(A,ES14.6,A)') " Total RMS (unweighted) = ", rms_total_fit * unit_con - write(std_out,'(A,ES14.6,A)') " Weighted RMS = ", rms_w_fit * unit_con - write(std_out,'(A,ES14.6,A)') " Energy RMS = ", en_rms_fit * unit_con - write(std_out,'(A,ES14.6,A)') " Dipole RMS = ", dip_rms_fit * unit_con - write(std_out, '(A)') " " - - write(std_out, '(A)') "VALIDATION / REFERENCE SET" - write(std_out, '(A)') "##################################" - write(std_out,'(A,ES14.6,A)') " Total RMS (unweighted) = ", rms_total_val * unit_con - write(std_out,'(A,ES14.6,A)') " Weighted RMS = ", rms_w_val * unit_con - write(std_out,'(A,ES14.6,A)') " Energy RMS = ", en_rms_val * unit_con - write(std_out,'(A,ES14.6,A)') " Dipole RMS = ", dip_rms_val * unit_con - write(std_out, '(A)') " " - - write(std_out, '(A)') "COMBINED (FIT + VALIDATION)" - write(std_out, '(A)') "------------------------------" - write(std_out,'(A,ES14.6,A)') " Total RMS (unweighted) = ", rms_total * unit_con - write(std_out,'(A,ES14.6,A)') " Weighted RMS = ", rms_weight * unit_con - write(std_out,'(A,ES14.6,A)') " Energy RMS = ", en_rms * unit_con - write(std_out,'(A,ES14.6,A)') " Dipole RMS = ", dip_rms * unit_con - write(std_out,'(A)')" END OF ERROR SUMMARY" - !write(std_out, '(A)') repeat("-", 70) - ! Optional: add more detailed per-component output if desired - ! do k = 1, inp_out - ! if (cnt(k) > 0) then - ! write(std_out,'(A,I0,A,ES12.5)') "Output ",k," RMS = ", & - ! sqrt(total_rms(k)/real(cnt(k),dp)) * unit_con - ! end if - ! end do - - end subroutine print_ErrorSummary - -end module print_error \ No newline at end of file