Merge pull request #11 from Trovemaster/hotfix/01

Hotfix/01

makefile

@@ -38,7 +38,7 @@ trove.x:        trove.o accuracy.o perturbation.o fields.o symmetry.o molecules.
 	            lapack.o plasma.o moltype.o refinement.o dipole.o refinement.o tran.o diag.o timer.o input.o \
                 mol_xy.o mol_xy2.o mol_xy3.o mol_xy4.o mol_zxy2.o mol_zxy3.o mol_ch3oh.o mol_abcd.o mol_c2h4.o mol_c2h6.o mol_c3h6.o \
 				pot_xy2.o pot_xy3.o pot_xy4.o pot_zxy2.o pot_zxy3.o pot_ch3oh.o pot_abcd.o pot_c2h4.o pot_c2h6.o pot_c3h6.o  $(pot_user).o \
-				prop_xy2.o prop_xy2_quad.o rotme_cart_tens.o fwigxjpf.o extfield.o richmol_data.o hyperfine.o
+				prop_xy2.o prop_xy2_quad.o prop_xy2_spinrot.o rotme_cart_tens.o fwigxjpf.o extfield.o richmol_data.o kin_xy2.o 
 				$(FOR) $(FFLAGS) -o j-trove$(PLAT).x $^ $(LIB)
 
 trove.o:        accuracy.o fields.o perturbation.o symmetry.o timer.o moltype.o dipole.o refinement.o tran.o extfield.o
@@ -47,7 +47,7 @@ fields.o:       accuracy.o molecules.o lapack.o me_str.o timer.o me_numer.o inpu
 symmetry.o:     accuracy.o
 molecules.o:    accuracy.o moltype.o mol_xy.o mol_xy2.o mol_xy3.o mol_xy4.o mol_zxy2.o mol_zxy3.o mol_ch3oh.o mol_abcd.o mol_c2h4.o mol_c3h6.o  mol_c2h6.o \
 				lapack.o pot_xy2.o pot_xy3.o mol_xy4.o pot_zxy2.o pot_zxy3.o pot_ch3oh.o pot_abcd.o pot_c2h4.o  pot_c3h6.o  pot_c2h6.o \
-				symmetry.o $(pot_user).o prop_xy2.o prop_xy2_quad.o
+				symmetry.o $(pot_user).o prop_xy2.o prop_xy2_quad.o prop_xy2_spinrot.o kin_xy2.o
 
 me_numer.o:     accuracy.o molecules.o timer.o
 me_str.o:       accuracy.o timer.o me_numer.o
@@ -97,10 +97,10 @@ kin_xy2.o:      accuracy.o moltype.o
 
 prop_xy2.o:     accuracy.o moltype.o timer.o pot_xy2.o
 prop_xy2_quad.o:     accuracy.o moltype.o timer.o pot_xy2.o
+prop_xy2_spinrot.o:  accuracy.o moltype.o timer.o pot_xy2.o
 
-extfield.o:         accuracy.o moltype.o fields.o symmetry.o tran.o timer.o rotme_cart_tens.o tran.o richmol_data.o hyperfine.o
+extfield.o:         accuracy.o moltype.o fields.o symmetry.o tran.o timer.o rotme_cart_tens.o tran.o richmol_data.o 
 rotme_cart_tens.o:  accuracy.o timer.o fwigxjpf.o
-hyperfine.o:        accuracy.o timer.o fwigxjpf.o richmol_data.o
 richmol_data.o:     accuracy.o timer.o
 fwigxjpf.o:
 	$(FOR) -c $(FFLAGS) -cpp fwigxjpf.f90

prop_xy2_quad.f90

@@ -1,3 +1,5 @@
+! Electric quadrupole moment tensor for XY2-type moleucles.
+
 module prop_xy2_quad
   use accuracy
   use moltype
@@ -19,6 +21,10 @@ module prop_xy2_quad
 !###################################################################################################################################
 
 
+! Electric quadrupole function in the bisector frame, here X axis bisects the
+! valence bond angle and the Z axis is aligned along the linear axis of the
+! molecule at linearity.
+
 recursive subroutine prop_xy2_qmom_bisect_frame(rank, ncoords, natoms, local, xyz, f)
 
   implicit none
@@ -83,7 +89,8 @@ end subroutine prop_xy2_qmom_bisect_frame
 !###################################################################################################################################
 
 
-! Subroutine to test xy2_qmom_bisect_frame if it is able to reproduce the original ab initio Cartesian components
+! This is to test the prop_xy2_qmom_bisect_frame functions if it is able to
+! reproduce the original data used for the fitting.
 
 subroutine TEST_xy2_qmom_bisect_frame(rank, ncoords, natoms, local, xyz, f)
 
@@ -393,4 +400,4 @@ function xy2_func_b2_n3_d6(coords, params) result(v)
 end function xy2_func_b2_n3_d6
 
 
-end module prop_xy2_quad
+end module prop_xy2_quad

xy2_spinrot.f90 → prop_xy2_spinrot.f90

@@ -1,4 +1,4 @@
-module xy2_spinrot
+module prop_xy2_spinrot
   use accuracy
   use moltype
   use pot_xy2, only : MLloc2pqr_xy2
@@ -23,6 +23,7 @@ subroutine prop_xy2_spin_rotation_bisector(rank, ncoords, natoms, local, xyz, f)
   real(ark) :: xyz0(3), xyz_(natoms,3),r1,r2,alpha,rho,rho_over_sinrho,rho2_over_sinrho2
   real(ark) :: c(3,3), mat(3,3), c_out(5,-2:0), e1(3), e2(3), e3(3),x(natoms,3)
   real(ark),parameter  :: rho_threshold = 0.01_rk
+  integer(ik) :: icentre
   !
   if (rank/=9) then
     write(out, '(/a,1x,i3,1x,a)') 'prop_xy2_spin_rotation_bisector: rank of the dipole moment vector =', rank, ', expected 9'
@@ -63,14 +64,27 @@ subroutine prop_xy2_spin_rotation_bisector(rank, ncoords, natoms, local, xyz, f)
   !
   rho = pi-alpha
   !
+  icentre = extF%coef(1,1) ! = 1 or 2, determines on which Y-atom spin-rotation tensor is centered
+  !
   ! fitted tensor elements
   c = 0
   !
-  c(1,1) = fit_xy2_sr(extF%nterms(1)-2, extF%coef(3:extF%nterms(1),1), (/r1, r2, alpha/))
-  c(2,2) = fit_xy2_sr(extF%nterms(2)-2, extF%coef(3:extF%nterms(2),2), (/r1, r2, alpha/))
-  c(3,3) = fit_xy2_sr(extF%nterms(3)-2, extF%coef(3:extF%nterms(3),3), (/r1, r2, alpha/))
-  c(1,3) = fit_xy2_sr(extF%nterms(4)-2, extF%coef(3:extF%nterms(4),4), (/r1, r2, alpha/))
-  c(3,1) = fit_xy2_sr(extF%nterms(5)-2, extF%coef(3:extF%nterms(5),5), (/r1, r2, alpha/))
+  if (icentre==1) then
+      c(1,1) = fit_xy2_sr(extF%nterms(1)-1, extF%coef(2:extF%nterms(1),1), (/r1, r2, alpha/))
+      c(2,2) = fit_xy2_sr(extF%nterms(2), extF%coef(1:extF%nterms(2),2), (/r1, r2, alpha/))
+      c(3,3) = fit_xy2_sr(extF%nterms(3), extF%coef(1:extF%nterms(3),3), (/r1, r2, alpha/))
+      c(1,3) = fit_xy2_sr(extF%nterms(4), extF%coef(1:extF%nterms(4),4), (/r1, r2, alpha/))
+      c(3,1) = fit_xy2_sr(extF%nterms(5), extF%coef(1:extF%nterms(5),5), (/r1, r2, alpha/))
+  elseif (icentre==2) then
+      c(1,1) = fit_xy2_sr(extF%nterms(1)-1, extF%coef(2:extF%nterms(1),1), (/r2, r1, alpha/))
+      c(2,2) = fit_xy2_sr(extF%nterms(2), extF%coef(1:extF%nterms(2),2), (/r2, r1, alpha/))
+      c(3,3) = fit_xy2_sr(extF%nterms(3), extF%coef(1:extF%nterms(3),3), (/r2, r1, alpha/))
+      c(1,3) = -fit_xy2_sr(extF%nterms(4), extF%coef(1:extF%nterms(4),4), (/r2, r1, alpha/))
+      c(3,1) = -fit_xy2_sr(extF%nterms(5), extF%coef(1:extF%nterms(5),5), (/r2, r1, alpha/))
+  else
+      write(out, '(a,1x,i3)') 'prop_xy2_qmom_bisect_frame error: icentre /= (1 or 2)'
+      stop 'prop_xy2_qmom_bisect_frame - bad icentre parameter' 
+  endif
   !
   ! inverse transform
   !
@@ -290,4 +304,4 @@ function xy2_func_n3_d6(coords, params) result(v)
 end function xy2_func_n3_d6
 
 
-end module xy2_spinrot
+end module prop_xy2_spinrot