new coordinate transformation for ND2H and NHD2 with dihedral

Trovemaster · Sep 7, 2023 · b8f7384 · b8f7384
1 parent c9d15f6
commit b8f7384
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 3 deletions.
diff --git a/mol_xy3.f90 b/mol_xy3.f90
@@ -334,11 +334,11 @@ subroutine ML_b0_XY3(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
                !
             enddo
             !
-            if (verbose>=6) then 
+            if (verbose>=3) then 
               write(out,"(i6)") molec%natoms
               !
               write(out,"(/'N',3x,3f14.8)") b0(1,:,i)
-              write(out,"( 'F',3x,3f14.8)") b0(2,:,i)
+              write(out,"( 'D',3x,3f14.8)") b0(2,:,i)
               write(out,"( 'H',3x,3f14.8)") b0(3,:,i)
               write(out,"( 'H',3x,3f14.8)") b0(4,:,i)
               !

diff --git a/mol_zxy2.f90 b/mol_zxy2.f90
@@ -449,7 +449,7 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
      real(ark),   intent(in),optional :: rho_borders(2)  ! rhomim, rhomax - borders
 
      real(ark)               :: a0(molec%Natoms,3),CM_shift,tau,cosalpha_2,costau,delta,cosrho,req(1:3),alphaeq(1:2),tau_
-     real(ark)               :: xieq(6),rho,theta,sint_2,theta12,beta
+     real(ark)               :: xieq(6),rho,theta,sint_2,theta12,beta,m1,m2,m3,m4,ax,ay
      integer(ik)             :: Nbonds,i,n
 
 
@@ -486,6 +486,11 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
          !
       end select 
       !
+      m1 = molec%AtomMasses(1) 
+      m2 = molec%AtomMasses(2) 
+      m3 = molec%AtomMasses(3) 
+      m4 = molec%AtomMasses(4)
+      !
       if (trim(molec%coords_transform)=='R-THETA-TAU-MEP') then 
          !
          !if (trim(molec%potentype)/='POTEN_ZXY2_MEP_R_ALPHA_RHO_POWERS') then
@@ -619,6 +624,20 @@ subroutine ML_b0_ZXY2(Npoints,Natoms,b0,rho_i,rho_ref,rho_borders)
                !
                !
             endif 
+            !
+            if ( abs(req(1)-req(2))<sqrt(small_).and.abs(req(1)-req(3))<sqrt(small_) ) then 
+               ax = sqrt( 1._ark-2.0_rk*cos(tau) )/(-1.0_ark+cos(tau))
+               ay= -cos(tau)/(-1.0_ark+cos(tau))
+               !
+               if (tau<0) ax = -ax
+               !
+               theta = atan2(ay,ax)
+               !
+               theta = acos( cos(tau)/(1.0_ark-cos(tau)) )
+               !
+            endif 
+
+
             !
             b0(1,1,i) = 0.0_ark
             b0(1,2,i) = 0.0_ark