diff --git a/src/model/invariants_no3.f90 b/src/model/invariants_no3.f90 deleted file mode 100644 index 2cfef00..0000000 --- a/src/model/invariants_no3.f90 +++ /dev/null @@ -1,82 +0,0 @@ -Module invariants_mod - implicit none - contains - !---------------------------------------------------- - subroutine invariants(a,xs,ys,xb,yb,b,inv) - implicit none - !include "nnparams.incl" - double precision, intent(in) :: a, xs, ys, xb, yb, b - double precision, intent(out) :: inv(4) - 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 - 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 )) - - ! the only non zero invariant for bend pure cuts - - inv(1) = invar(1) - inv(2) = invar(5) - inv(3) = invar(9) - inv(4) = invar(23) - 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 -end module invariants_mod diff --git a/src/model/model_nh3.f90 b/src/model/model_nh3.f90 deleted file mode 100644 index b261634..0000000 --- a/src/model/model_nh3.f90 +++ /dev/null @@ -1,401 +0,0 @@ -module diabmodel - use iso_fortran_env, only: idp => int32, dp => real64 - use dip_param - implicit none - include "nnparams.incl" - integer(idp),parameter:: ndiab=4 - logical :: debug=.false. - 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(ndiab,ndiab) - integer(idp) id,i,j ,ii - real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) - real(dp),dimension(12):: shift,scal - xs=q(6) - ys=q(7) - xb=q(8) - yb=q(9) - a=q(5) - b=q(10) - - call init_dip_planar_data() - - ! modify the parametr - !shift=1.0_dp - !scal=1.0d-2 - scal = [0.69753,0.44797,51.14259,2.76924,1.45300,91.58246, & - 14.07390,1.02550,1.68623,4.80804,7.69958,0.97871] - shift = [0.23479,0.26845,35.05940,2.27175,-0.33017,117.48895, & - -1.68211,0.79418,-1.60443,-10.41309,8.47695,1.25334] - ! V term of A2'' - ii=1 - do i =1,np - if (p(i) .ne. 0.0d0) then - p(i) =p(i)*(shift(ii) + scal(ii)*nn_out(ii) ) - ii=ii+1 - else - p(i) = p(i) - endif - enddo - - 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=1 !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(ndiab,ndiab) - integer(idp) id,i,j, ii - real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8) - real(dp),dimension(12):: shift,scal - xs=q(6) - ys=q(7) - xb=q(8) - yb=q(9) - a=q(5) - b=q(10) - - call init_dip_planar_data() - ! modify the parametr - !shift=1.0_dp - !scal=1.0d-3 - scal = [0.69753,0.44797,51.14259,2.76924,1.45300,91.58246, & - 14.07390,1.02550,1.68623,4.80804,7.69958,0.97871] - shift = [0.23479,0.26845,35.05940,2.27175,-0.33017,117.48895, & - -1.68211,0.79418,-1.60443,-10.41309,8.47695,1.25334] - - ! V term of A2'' - ii=1 - do i =1,np - if (p(i) .ne. 0) then - p(i) =p(i)*(shift(ii) + scal(ii)*nn_out(ii) ) - ii=ii+1 - else - p(i) = p(i) - endif - enddo - - - - 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=1 !1 - e(1,1)=e(1,1)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb - - id=id+1 !2 - e(2,2)=e(2,2)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb - e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb - id=id+1 !3 - e(4,4)=e(4,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb - ! order 2 - id=id+1 !4 - e(1,1)=e(1,1)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & - -p(pst(1,id)+2)*(xs*yb+xb*ys) - id=id+1 !5 - e(2,2)=e(2,2)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & - -p(pst(1,id)+2)*(xs*yb+xb*ys) - - - e(3,3)=e(3,3)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & - -p(pst(1,id)+2)*(xs*yb+xb*ys) - id=id+1 !6 - e(4,4)=e(4,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) & - -p(pst(1,id)+2)*(xs*yb+xb*ys) - ! order 3 - id=id+1 !7 - e(1,1)=e(1,1)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb - id=id+1 !8 - e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb - e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb - id=id+1 !9 - e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb - - ! V- term + totally symmetric coord a - - ! JAHN TELLER COUPLING TERM - ! order 0 - id=id+1 !10 - e(2,3)=e(2,3)+p(pst(1,id)) - ! order 1 - - id=id+1 !11 - e(2,2)=e(2,2)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb - e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb - e(2,3)=e(2,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb - !id=id+1 !12 - ! order 2 - id=id+1 !12 - e(2,2)=e(2,2)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+p(pst(1,id)+2)*(xs*yb+xb*ys) - e(3,3)=e(3,3)-p(pst(1,id))*2*xs*ys-p(pst(1,id)+1)*2*xb*yb-p(pst(1,id)+2)*(xs*yb+xb*ys) - e(2,3)=e(2,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)) & - -p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb - ! order 3 - id=id+1 !13 - do i=1,4 - j=i-1 - e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) - e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4) - e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i) - enddo - e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb - e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb - e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb - ! PSEUDO JAHN TELLER - ! ORDER 0 - ! THE COUPLING OF A2 GROUND STATE WITH E - ! ################################################### - ! ################################################### - ! order 0 - id=id+1 !14 - e(1,3)=e(1,3)-b*(p(pst(1,id))) - ! order 1 - id=id+1 !15 - e(1,2)=e(1,2)-b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) - e(1,3)=e(1,3)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb) - - ! order 2 - id=id+1 !16 - e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& - +p(pst(1,id)+2)*(xs*yb+xb*ys)) - e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)& - +p(pst(1,id)+2)*(xs*xb-ys*yb) - p(pst(1,id)+3)*(xs**2+ys**2)) - - ! order 3 - id = id+1 ! 17 - do i=1,4 - e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) - e(1,3)=e(1,3)-b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) - enddo - - - ! THE COUPLING OF A2 WITH A1 - !#################################################### - !#################################################### - ! order 1 - id=id+1 !17 - e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb) - ! order 2 - id=id+1 !18 - e(1,4)=e(1,4)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)& - +p(pst(1,id)+2)*(xs*yb+xb*ys)) - - - ! THE COUPLING OF A1 WITH E - !#################################################### - !#################################################### - ! order 0 - id=id+1 !19 - e(3,4)=e(3,4)-p(pst(1,id)) - ! order 1 - id=id+1 !20 - e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb - e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb - ! order 2 - id=id+1 !21 - e(2,4)=e(2,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb) & - +p(pst(1,id)+2)*(xs*yb+xb*ys) - e(3,4)=e(3,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) & - +p(pst(1,id)+2)*(xs*xb-ys*yb) - p(pst(1,id)+3)*(xs**2+ys**2) - - id =id+1 ! 23 - ! order 3 - do i=1,4 - e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4) - e(3,4)=e(3,4)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i) - enddo - ! end of the model - e(2,1)=e(1,2) - e(3,1)=e(1,3) - e(3,2)=e(2,3) - e(4,1)=e(1,4) - e(4,2)=e(2,4) - e(4,3)=e(3,4) -end subroutine diab_y - - subroutine copy_2_lower_triangle(mat) - real(dp), intent(inout) :: mat(:, :) - integer :: m, n -! write lower triangle of matrix symmetrical - do n = 2, size(mat, 1) - do m = 1, n - 1 - mat(n, m) = mat(m, n) - end do - end do - end subroutine copy_2_lower_triangle - -end module diabmodel diff --git a/src/model/model_nh3_py.f90 b/src/model/model_nh3_py.f90 index a919f96..ebe7b1a 100644 --- a/src/model/model_nh3_py.f90 +++ b/src/model/model_nh3_py.f90 @@ -24,12 +24,14 @@ module diabmodel call init_dip_planar_data() ! modify the parametr - shift=1.0_dp - scal=1.0d-2 - !scal = [0.69753,0.44797,51.14259,2.76924,1.45300,91.58246, & - ! 14.07390,1.02550,1.68623,4.80804,7.69958,0.97871] - !shift = [0.23479,0.26845,35.05940,2.27175,-0.33017,117.48895, & - ! -1.68211,0.79418,-1.60443,-10.41309,8.47695,1.25334] + !shift=1.0_dp + !scal=1.0d-2 + scal = [81.50814,1.19850,245.66937,885.73515,875.91907,4.55261,4.14318,17.98337, & + 47.60657,24.68377,124.20271,40.13721,153.81484,13.19697,258.11015,769.62299, & + 147.59778,97.94602,413.29279,62.33428,1.88373] + shift = [-45.54845,-0.12528,139.98528,887.54027,602.90180,-1.61771,2.28343,-12.43082, & + -29.81680,9.29555,123.83516,-16.80044,108.07947,-2.19734,120.51791,496.11692, & + -83.20774,43.09530,373.00489,-15.32774,1.30491] ! V term of A2'' ii=1 do i =1,np @@ -216,13 +218,15 @@ module diabmodel call init_dip_planar_data() ! modify the parametr - shift=1.0_dp - scal=1.0d-2 - !scal = [0.69753,0.44797,51.14259,2.76924,1.45300,91.58246, & - ! 14.07390,1.02550,1.68623,4.80804,7.69958,0.97871] - !shift = [0.23479,0.26845,35.05940,2.27175,-0.33017,117.48895, & - ! -1.68211,0.79418,-1.60443,-10.41309,8.47695,1.25334] - + !shift=1.0_dp + !scal=1.0d-2 + scal = [81.50814,1.19850,245.66937,885.73515,875.91907,4.55261,4.14318,17.98337, & + 47.60657,24.68377,124.20271,40.13721,153.81484,13.19697,258.11015,769.62299, & + 147.59778,97.94602,413.29279,62.33428,1.88373] + shift = [-45.54845,-0.12528,139.98528,887.54027,602.90180,-1.61771,2.28343,-12.43082, & + -29.81680,9.29555,123.83516,-16.80044,108.07947,-2.19734,120.51791,496.11692, & + -83.20774,43.09530,373.00489,-15.32774,1.30491] + ! V term of A2'' ii=1 do i =1,np diff --git a/src/model/nncoords_no3.f90 b/src/model/nncoords_no3.f90 deleted file mode 100644 index 0005ee0..0000000 --- a/src/model/nncoords_no3.f90 +++ /dev/null @@ -1,716 +0,0 @@ - -!-------------------------------------------------------------------------------------- - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % SUBROUTINE CTRANS(...) -! % -! % M. Vossel 21.03.2023 -! % -! % Routine to transform symmetryinput coordinates to symmetrized -! % coordinates. Distances Are discribet by Morse coordinates or -! % TMC depending on Set Parameters in the Genetic Input. -! % -! % input variables -! % q: -! % q(1): H1x -! % q(2): y -! % q(3): z -! % q(4): H2x -! % q(5): y -! % q(6): z -! % q(7): H3x -! % q(8): y -! % q(9): z -! -! -! -! % Internal variables: -! % t: primitive coordinates (double[qn]) -! % t(1): -! % t(2): -! % t(3): -! % t(4): -! % t(5): -! % t(6): -! % t(7): -! % t(8): -! % t(9): -! % t: dummy (double[qn]) -! % p: parameter vector -! % npar: length of parameter vector -! % -! % Output variables -! % s: symmetrized coordinates (double[qn]) -! % s(1): CH-symetric streatch -! % s(2): CH-asymetric streatch-ex -! % s(3): CH-asymetric streatch-ey -! % s(4): CH-bend-ex -! % s(5): CH-bend-ey -! % s(6): CH-umbrella -! % s(7): CH-umbrella**2 -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -module ctrans_mod - implicit none - !include 'only_model.incl' - include 'nnparams.incl' -! precalculate pi, 2*pi and angle to radian conversion - double precision, parameter :: pii = 4.00d0*datan(1.00d0) - double precision, parameter :: pi2 = 2.00d0*pii - double precision, parameter :: ang2rad = pii/180.00d0 -! precalculate roots - double precision, parameter:: sq2 = 1.00d0/dsqrt(2.00d0) - double precision, parameter:: sq3 = 1.00d0/dsqrt(3.00d0) - double precision, parameter:: sq6 = 1.00d0/dsqrt(6.00d0) -! change distances for equilibrium - double precision, parameter :: dchequi = 1.0228710942d0 ! req NH3+ - -contains - subroutine ctrans(q) - !use dim_parameter, only: qn - use invariants_mod, only: invariants - integer k !running indices - double precision, intent(inout) :: q(maxnin) !given coordinates - double precision :: s(maxnin) !output coordinates symmetry adapted and scaled - double precision :: t(maxnin) !output coordinates symmetry adapted but not scaled - ! ANN Variables - !double precision, optional, intent(out) :: invariants(:) - ! kartesian coordianates copy from MeF+ so substitute c by n and removed f - double precision ch1(3), ch2(3), ch3(3), c_atom(3) - double precision nh1(3), nh2(3), nh3(3) - double precision zaxis(3), xaxis(3), yaxis(3) - double precision ph1(3), ph2(3), ph3(3) - ! primitive coordinates - double precision dch1, dch2, dch3 !nh-distances - double precision umb !Umbrella Angle from xy-plane - - ! Symmetry coordinates - double precision aR !a1-modes H-Dist., - double precision exR, exAng !ex components H-Dist., H-Ang. - double precision eyR, eyAng !ey components H-Dist., H-Ang. - double precision inv(4) - ! debugging - logical, parameter :: dbg = .false. - - ! initialize coordinate vectors - s = 0.0d0 - t = 0.0d0 - - ! write kartesian coords for readability - c_atom(1:3) = 0.0d0 ! N-atom at origin - do k = 1, 3 - ch1(k) = q(k ) - ch2(k) = q(k + 3) - ch3(k) = q(k + 6) - end do - q=0.d0 - - ! construct z-axis - nh1 = normalized(ch1) - nh2 = normalized(ch2) - nh3 = normalized(ch3) - zaxis = create_plane(nh1, nh2, nh3) - - ! calculate bonding distance - dch1 = norm(ch1) - dch2 = norm(ch2) - dch3 = norm(ch3) - - ! construct symmertic and antisymmetric strech - aR = symmetrize(dch1 - dchequi, dch2 - dchequi, dch3 - dchequi, 'a') - exR = symmetrize(dch1, dch2, dch3, 'x') - eyR = symmetrize(dch1, dch2, dch3, 'y') - - ! construc x-axis and y axis - ph1 = normalized(project_point_into_plane(nh1, zaxis, c_atom)) - xaxis = normalized(ph1) - yaxis = xproduct(zaxis, xaxis) ! right hand side koordinates - - ! project H atoms into C plane - ph2 = normalized(project_point_into_plane(nh2, zaxis, c_atom)) - ph3 = normalized(project_point_into_plane(nh3, zaxis, c_atom)) - - call construct_HBend(exAng, eyAng, ph1, ph2, ph3, xaxis, yaxis) - umb = construct_umbrella(nh1, nh2, nh3, zaxis) - - ! set symmetry coordinates and even powers of umbrella - !s(1) = dch1-dchequi!aR - !s(2) = dch2-dchequi!exR - !s(3) = dch3-dchequi!eyR - ! call invariants and get them - ! 24 invariants - - call invariants(aR,exR,eyR,exAng,eyAng,umb,inv) - q(1:4)=inv(1:4) - q(5) = aR - q(6) = exR - q(7) = eyR - q(8) = exAng - q(9) = eyAng - q(10) = umb - ! pairwise distances as second coordinate set - !call pair_distance(q, t(1:6)) - - !if (dbg) write (6, '("sym coords s=",9f16.8)') s(1:qn) - !if (dbg) write (6, '("sym coords t=",9f16.8)') t(1:qn) - !if (present(invariants)) then - ! call get_invariants(s, invariants) - !end if - ! RETURN Q AS INTERNAL COORD - end subroutine ctrans -! subroutine ctrans(q) -! use invariants_mod, only: invariants -! implicit none -! double precision, intent(inout):: q(maxnin) -! double precision:: invar(24) -! double precision:: a,b,esx,esy,ebx,eby -! a=q(1) -! esx=q(2) -! esy=q(3) -! ebx=q(4) -! eby=q(5) -! b=q(6) -! call invariants(a,esx,esy,ebx,eby,b,invar) -! -! q(1:24)=invar(1:24) -! q(25)=esx -! q(26)=esy -! q(27)=ebx -! q(28)=eby -! q(29)=b -! end subroutine ctrans - subroutine pair_distance(q, r) - double precision, intent(in) :: q(9) - double precision, intent(out) :: r(6) - double precision :: atom(3, 4) - integer :: n, k, count - - !atom order: H1 H2 H3 N - atom(:, 1:3) = reshape(q, [3, 3]) - atom(:, 4) = (/0.00d0, 0.00d0, 0.00d0/) - - ! disntace order 12 13 14 23 24 34 - count = 0 - do n = 1, size(atom, 2) - do k = n + 1, size(atom, 2) - count = count + 1 - r(count) = sqrt(sum((atom(:, k) - atom(:, n))**2)) - end do - end do - end subroutine pair_distance - - function morse_and_symmetrize(x,p,pst) result(s) - double precision, intent(in),dimension(3) :: x - double precision, intent(in),dimension(11) :: p - integer, intent(in),dimension(2) :: pst - integer :: k - double precision, dimension(3) :: s - double precision, dimension(3) :: t - - ! Morse transform - do k=1,3 - t(k) = morse_transform(x(k), p, pst) - end do - s(1) = symmetrize(t(1), t(2), t(3), 'a') - s(2) = symmetrize(t(1), t(2), t(3), 'x') - s(3) = symmetrize(t(1), t(2), t(3), 'y') - end function morse_and_symmetrize - -! subroutine get_invariants(s, inv_out) -! use dim_parameter, only: qn -! use select_monom_mod, only: v_e_monom, v_ee_monom -! double precision, intent(in) :: s(qn) -! double precision, intent(out) :: inv_out(:) -! ! double precision, parameter :: ck = 1.00d0, dk = 1.00d0/ck ! scaling for higher order invariants -! double precision inv(24) -! integer, parameter :: inv_order(12) = & ! the order in which the invariants are selected -! & [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] -! double precision Rch, umb, xR, yR, xAng, yAng -!! for readability -! Rch = s(1) -! xR = s(2) -! yR = s(3) -! xAng = s(4) -! yAng = s(5) -! umb = s(6)**2 -!! invarianten -!! a moden -! inv(1) = Rch -! inv(2) = umb -!! invariante e pairs -! inv(3) = v_e_monom(xR, yR, 1) -! inv(4) = v_e_monom(xAng, yAng, 1) -!! third order e pairs -! inv(5) = v_e_monom(xR, yR, 2) -! inv(6) = v_e_monom(xAng, yAng, 2) -! ! invariant ee coupling -! inv(7) = v_ee_monom(xR, yR, xAng, yAng, 1) -! ! mode combinations -! inv(8) = Rch*umb -! -! inv(9) = Rch*v_e_monom(xR, yR, 1) -! inv(10) = umb*v_e_monom(xR, yR, 1) -! -! inv(11) = Rch*v_e_monom(xAng, yAng, 1) -! inv(12) = umb*v_e_monom(xAng, yAng, 1) -! -!! damp coordinates because of second order and higher invariants -! inv(3) = sign(sqrt(abs(inv(3))), inv(3)) -! inv(4) = sign(sqrt(abs(inv(4))), inv(4)) -! inv(5) = sign((abs(inv(5))**(1./3.)), inv(5)) -! inv(6) = sign((abs(inv(6))**(1./3.)), inv(6)) -! inv(7) = sign((abs(inv(7))**(1./3.)), inv(7)) -! inv(8) = sign(sqrt(abs(inv(8))), inv(8)) -! inv(9) = sign((abs(inv(9))**(1./3.)), inv(9)) -! inv(10) = sign((abs(inv(10))**(1./3.)), inv(10)) -! inv(11) = sign((abs(inv(11))**(1./3.)), inv(11)) -! inv(12) = sign((abs(inv(12))**(1./3.)), inv(12)) -! -! inv_out(:) = inv(inv_order(1:size(inv_out, 1))) -! -! end subroutine get_invariants - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % real part of spherical harmonics -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! Ylm shifted to 0 for theta=0 -! double precision function ylm(theta, phi, l, m) -! implicit none -! double precision theta, phi -! integer l, m -! ylm = plm2(dcos(theta), l, m)*cos(m*phi) - plm2(1.00d0, l, m) -! end function ylm -!!---------------------------------------------------------- -! double precision function plm2(x, l, n) -! implicit none -! double precision x -! integer l, m, n -! -! double precision pmm, p_mp1m, pllm -! integer ll -! -!! negative m und bereich von x abfangen -! if ((l .lt. 0)& -! &.or. (abs(n) .gt. abs(l))& -! &.or. (abs(x) .gt. 1.)) then -! write (6, '(''bad arguments in legendre'')') -! stop -! end if -! -!! fix sign of m to compute the positiv m -! m = abs(n) -! -! pmm = (-1)**m*dsqrt(fac(2*m))*1./((2**m)*fac(m))& !compute P(m,m) not P(l,l) -! &*(dsqrt(1.-x**2))**m -! -! if (l .eq. m) then -! plm2 = pmm !P(l,m)=P(m,m) -! else -! p_mp1m = x*dsqrt(dble(2*m + 1))*pmm !compute P(m+1,m) -! if (l .eq. m + 1) then -! plm2 = p_mp1m !P(l,m)=P(m+1,m) -! else -! do ll = m + 2, l -! pllm = x*(2*ll - 1)/dsqrt(dble(ll**2 - m**2))*p_mp1m& ! compute P(m+2,m) up to P(l,m) recursively -! &- dsqrt(dble((ll - 1)**2 - m**2))& -! &/dsqrt(dble(l**2 - m**2))*pmm -!! schreibe m+2 und m+1 jeweils fuer die naechste iteration -! pmm = p_mp1m !P(m,m) = P(m+1,m) -! p_mp1m = pllm !P(m+1,m) = P(m+2,m) -! end do -! plm2 = pllm !P(l,m)=P(m+k,m), k element N -! end if -! end if -! -!! sets the phase of -m term right (ignored to gurantee Ylm=(Yl-m)* for JT terms -!! if(n.lt.0) then -!! plm2 = (-1)**m * plm2 !* fac(l-m)/fac(l+m) -!! endif -! -! end function -!---------------------------------------------------------------------------------------------------- - double precision function fac(i) - integer i - select case (i) - case (0) - fac = 1.00d0 - case (1) - fac = 1.00d0 - case (2) - fac = 2.00d0 - case (3) - fac = 6.00d0 - case (4) - fac = 24.00d0 - case (5) - fac = 120.00d0 - case (6) - fac = 720.00d0 - case (7) - fac = 5040.00d0 - case (8) - fac = 40320.00d0 - case (9) - fac = 362880.00d0 - case (10) - fac = 3628800.00d0 - case (11) - fac = 39916800.00d0 - case (12) - fac = 479001600.00d0 - case default - write (*, *) 'ERROR: no case for given faculty, Max is 12!' - stop - end select - end function fac - - ! Does the simplest morse transform possible - ! one skaling factor + shift - function morse_transform(x, p, pst) result(t) - double precision, intent(in) :: x - double precision, intent(in) :: p(11) - integer, intent(in) :: pst(2) - double precision :: t - if (pst(2) == 11) then - t = 1.00d0 - exp(-abs(p(2))*(x - p(1))) - else - error stop 'in morse_transform key required or wrong number of parameters' - end if - end function morse_transform - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % FUNCTION F(...) ! MAIK DEPRICATING OVER THE TOP MORSE FUNCTION FOR MYSELF -! % -! % Returns exponent of tunable Morse coordinate -! % exponent is polynomial * gaussian (skewed) -! % ilabel = 1 or 2 selects the parameters a and sfac to be used -! % -! % Background: better representation of the prefector in the -! % exponend of the morse function. -! % Formular: f(r) = lest no3 paper -! % -! % Variables: -! % x: distance of atoms (double) -! % p: parameter vector (double[20]) -! % ii: 1 for CCl and 2 for CCH (int) -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - pure function f(x, p, ii) - integer , intent(in) :: ii !1 for CCL and 2 for CCH - double precision, intent(in) :: x !coordinate - double precision, intent(in) :: p(11) !parameter-vector - - integer i !running index - - double precision r !equilibrium distance - double precision gaus !gaus part of f - double precision poly !polynom part of f - double precision skew !tanh part of f - - double precision f !prefactor of exponent and returned value - - integer npoly(2) !order of polynom - -! Maximum polynom order - npoly(1) = 5 - npoly(2) = 5 - -! p(1): position of equilibrium -! p(2): constant of exponent -! p(3): constant for skewing the gaussian -! p(4): tuning for skewing the gaussian -! p(5): Gaussian exponent -! p(6): Shift of Gaussian maximum -! p(7)...: polynomial coefficients -! p(8+n)...: coefficients of Morse Power series - -! 1-exp{[p(2)+exp{-p(5)[x-p(6)]^2}[Taylor{p(7+n)}(x-p(6))]][x-p(1)]} - -! Tunable Morse function -! Power series in Tunable Morse coordinates of order m -! exponent is polynomial of order npoly * gaussian + switching function - -! set r r-r_e - r = x - r = r - p(1) - -! set up skewing function: - skew = 0.50d0*p(3)*(dtanh(dabs(p(4))*(r - p(6))) + 1.00d0) - -! set up gaussian function: - gaus = dexp(-dabs(p(5))*(r - p(6))**2) - -! set up power series: - poly = 0.00d0 - do i = 0, npoly(ii) - 1 - poly = poly + p(7 + i)*(r - p(6))**i - end do -! set up full exponent function: - f = dabs(p(2)) + skew + gaus*poly - - end function -!---------------------------------------------------------------------------------------------------- - pure function xproduct(a, b) result(axb) - double precision, intent(in) :: a(3), b(3) - double precision :: axb(3) !crossproduct a x b - axb(1) = a(2)*b(3) - a(3)*b(2) - axb(2) = a(3)*b(1) - a(1)*b(3) - axb(3) = a(1)*b(2) - a(2)*b(1) - end function xproduct - - pure function normalized(v) result(r) - double precision, intent(in) :: v(:) - double precision :: r(size(v)) - r = v/norm(v) - end function normalized - - pure function norm(v) result(n) - double precision, intent(in) :: v(:) - double precision n - n = dsqrt(sum(v(:)**2)) - end function norm - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % FUNCTION Project_Point_Into_Plane(x,n,r0) result(p) -! % return the to n orthogonal part of a vector x-r0 -! % p: projected point in plane -! % x: point being projected -! % n: normalvector of plane -! % r0: Point in plane -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - pure function project_point_into_plane(x, n, r0) result(p) - double precision, intent(in) :: x(:), n(:), r0(:) - double precision :: p(size(x)), xs(size(x)) - xs = x - r0 - p = xs - plane_to_point(x, n, r0) - end function project_point_into_plane - -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -! % Function Plane_To_Point(x,n,r0) result(p) -! % p: part of n in x -! % x: point being projected -! % n: normalvector of plane -! % r0: Point in plane -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - pure function plane_to_point(x, n, r0) result(p) - double precision, intent(in) :: x(:), n(:), r0(:) - double precision p(size(x)), xs(size(x)), nn(size(n)) - nn = normalized(n) - xs = x - r0 - p = dot_product(nn, xs)*nn - end function plane_to_point - - subroutine check_coordinates(q) -! check for faulty kartesain coordinates - double precision, intent(in) :: q(:) - integer :: i - if (all(abs(q) <= epsilon(0.00d0))) then - stop 'Error (ctrans): all kartesian coordinates are<=1d-8' - end if - do i = 1, 9, 3 - if (all(abs(q(i:i + 2)) <= epsilon(0.00d0))) then - write (*, *) q - stop 'Error(ctrans):kartesian coordinates zero for one atom' - end if - end do - end subroutine - - pure function rotor_a_to_z(a, z) result(r) - double precision, intent(in) :: a(3), z(3) - double precision :: r(3, 3) - double precision :: alpha - double precision :: s1(3), s(3, 3), rotor(3, 3) - s1 = xproduct(normalized(a), normalized(z)) - alpha = asin(norm(s1)) - s(:, 1) = normalized(s1) - s(:, 2) = normalized(z) - s(:, 3) = xproduct(s1, z) - rotor = init_rotor(alpha, 0.00d0, 0.00d0) - r = matmul(s, matmul(rotor, transpose(s))) - end function - -! function returning Rz(gamma) * Ry(beta) * Rx(alpha) for basis order xyz - pure function init_rotor(alpha, beta, gamma) result(rotor) - double precision, intent(in) :: alpha, beta, gamma - double precision :: rotor(3, 3) - rotor = 0.00d0 - rotor(1, 1) = dcos(beta)*dcos(gamma) - rotor(1, 2) = dsin(alpha)*dsin(beta)*dcos(gamma)& - &- dcos(alpha)*dsin(gamma) - rotor(1, 3) = dcos(alpha)*dsin(beta)*dcos(gamma)& - &+ dsin(alpha)*dsin(gamma) - - rotor(2, 1) = dcos(beta)*dsin(gamma) - rotor(2, 2) = dsin(alpha)*dsin(beta)*dsin(gamma)& - &+ dcos(alpha)*dcos(gamma) - rotor(2, 3) = dcos(alpha)*dsin(beta)*dsin(gamma)& - &- dsin(alpha)*dcos(gamma) - - rotor(3, 1) = -dsin(beta) - rotor(3, 2) = dsin(alpha)*dcos(beta) - rotor(3, 3) = dcos(alpha)*dcos(beta) - end function init_rotor - - pure function create_plane(a, b, c) result(n) - double precision, intent(in) :: a(3), b(3), c(3) - double precision :: n(3) - double precision :: axb(3), bxc(3), cxa(3) - axb = xproduct(a, b) - bxc = xproduct(b, c) - cxa = xproduct(c, a) - n = normalized(axb + bxc + cxa) - end function create_plane - - function symmetrize(q1, q2, q3, sym) result(s) - double precision, intent(in) :: q1, q2, q3 - character, intent(in) :: sym - double precision :: s - select case (sym) - case ('a') - s = (q1 + q2 + q3)*sq3 - case ('x') - s = sq6*(2.00d0*q1 - q2 - q3) - case ('y') - s = sq2*(q2 - q3) - case default - write (*, *) 'ERROR: no rule for symmetrize with sym=', sym - stop - end select - end function symmetrize - - subroutine construct_HBend(ex, ey, ph1, ph2, ph3, x_axis, y_axis) - double precision, intent(in) :: ph1(3), ph2(3), ph3(3) - double precision, intent(in) :: x_axis(3), y_axis(3) - double precision, intent(out) :: ex, ey - double precision :: x1, y1, alpha1 - double precision :: x2, y2, alpha2 - double precision :: x3, y3, alpha3 -! get x and y components of projected points - x1 = dot_product(ph1, x_axis) - y1 = dot_product(ph1, y_axis) - x2 = dot_product(ph2, x_axis) - y2 = dot_product(ph2, y_axis) - x3 = dot_product(ph3, x_axis) - y3 = dot_product(ph3, y_axis) -! -> calculate H deformation angles - alpha3 = datan2(y2, x2) - alpha2 = -datan2(y3, x3) !-120*ang2rad -! write(*,*)' atan2' -! write(*,*) 'alpha2:' , alpha2/ang2rad -! write(*,*) 'alpha3:' , alpha3/ang2rad - if (alpha2 .lt. 0) alpha2 = alpha2 + pi2 - if (alpha3 .lt. 0) alpha3 = alpha3 + pi2 - alpha1 = (pi2 - alpha2 - alpha3) -! write(*,*)' fixed break line' -! write(*,*) 'alpha1:' , alpha1/ang2rad -! write(*,*) 'alpha2:' , alpha2/ang2rad -! write(*,*) 'alpha3:' , alpha3/ang2rad - alpha1 = alpha1 !- 120.00d0*ang2rad - alpha2 = alpha2 !- 120.00d0*ang2rad - alpha3 = alpha3 !- 120.00d0*ang2rad -! write(*,*)' delta alpha' -! write(*,*) 'alpha1:' , alpha1/ang2rad -! write(*,*) 'alpha2:' , alpha2/ang2rad -! write(*,*) 'alpha3:' , alpha3/ang2rad -! write(*,*) - -! construct symmetric and antisymmetric H angles - ex = symmetrize(alpha1, alpha2, alpha3, 'x') - ey = symmetrize(alpha1, alpha2, alpha3, 'y') - end subroutine construct_HBend - - pure function construct_umbrella(nh1, nh2, nh3, n)& - &result(umb) - double precision, intent(in) :: nh1(3), nh2(3), nh3(3) - double precision, intent(in) :: n(3) - double precision :: umb - double precision :: theta(3) -! calculate projections for umberella angle - theta(1) = dacos(dot_product(n, nh1)) - theta(2) = dacos(dot_product(n, nh2)) - theta(3) = dacos(dot_product(n, nh3)) -! construct umberella angle - umb = sum(theta(1:3))/3.00d0 - 90.00d0*ang2rad - end function construct_umbrella - - pure subroutine construct_sphericals& - &(theta, phi, cf, xaxis, yaxis, zaxis) - double precision, intent(in) :: cf(3), xaxis(3), yaxis(3), zaxis(3) - double precision, intent(out) :: theta, phi - double precision :: x, y, z, v(3) - v = normalized(cf) - x = dot_product(v, normalized(xaxis)) - y = dot_product(v, normalized(yaxis)) - z = dot_product(v, normalized(zaxis)) - theta = dacos(z) - phi = -datan2(y, x) - end subroutine construct_sphericals - -! subroutine int2kart(internal, kart) -! double precision, intent(in) :: internal(6) -! double precision, intent(out) :: kart(9) -! double precision :: h1x, h1y, h1z -! double precision :: h2x, h2y, h2z -! double precision :: h3x, h3y, h3z -! double precision :: dch0, dch1, dch2, dch3 -! double precision :: a1, a2, a3, wci -! -! kart = 0.00d0 -! dch1 = dchequi + sq3*internal(1) + 2*sq6*internal(2) -! dch2 = dchequi + sq3*internal(1) - sq6*internal(2) + sq2*internal(3) -! dch3 = dchequi + sq3*internal(1) - sq6*internal(2) - sq2*internal(3) -! a1 = 2*sq6*internal(4) -! a2 = -sq6*internal(4) + sq2*internal(5) -! a3 = -sq6*internal(4) - sq2*internal(5) -! wci = internal(6) -! -! ! Berechnung kartesische Koordinaten -! ! ----------------------- -! h1x = dch1*cos(wci*ang2rad) -! h1y = 0.0 -! h1z = -dch1*sin(wci*ang2rad) -! -! h3x = dch2*cos((a2 + 120)*ang2rad)*cos(wci*ang2rad) -! h3y = dch2*sin((a2 + 120)*ang2rad)*cos(wci*ang2rad) -! h3z = -dch2*sin(wci*ang2rad) -! -! h2x = dch3*cos((-a3 - 120)*ang2rad)*cos(wci*ang2rad) -! h2y = dch3*sin((-a3 - 120)*ang2rad)*cos(wci*ang2rad) -! h2z = -dch3*sin(wci*ang2rad) -! -! kart(1) = h1x -! kart(2) = h1y -! kart(3) = h1z -! kart(4) = h2x -! kart(5) = h2y -! kart(6) = h2z -! kart(7) = h3x -! kart(8) = h3y -! kart(9) = h3z -! end subroutine int2kart - -end module ctrans_mod -!**** Define coordinate transformation applied to the input before fit. -!*** -!*** -!****Conventions: -!*** -!*** ctrans: subroutine transforming a single point in coordinate space - - subroutine trans_in(pat_in,ntot) - use ctrans_mod, only: ctrans - implicit none - - include 'nnparams.incl' - !include 'nndbg.incl' - !include 'nncommon.incl' - - double precision pat_in(maxnin,maxpats) - integer ntot - - integer j - - do j=1,ntot - call ctrans(pat_in(:,j)) - ! FIRST ELEMENT OF PAT-IN ARE USED BY NEURON NETWORK - write(62,'(6f16.8)') pat_in(4:9,j) - enddo - end