src/alphas.f90: corrections
authorAnthony Stone <ajs1@cam.ac.uk>
Tue, 16 Oct 2018 15:25:23 +0000 (16:25 +0100)
committerAnthony Stone <ajs1@cam.ac.uk>
Tue, 16 Oct 2018 15:25:23 +0000 (16:25 +0100)
read_site_pairs didn't allow for in-line polarizability input.
total_dd now prints mean polarizability and anisotropy as well as the
full polarizability tensor.

src/alphas.f90

index c1e40f6..741ef23 100644 (file)
@@ -1012,6 +1012,8 @@ END SUBROUTINE read_alphas
 
 SUBROUTINE read_site_pairs(lmax)
 
+!  Read polarizabilities in new format
+
 INTEGER, INTENT(OUT) :: lmax
 
 CHARACTER(LEN=20) :: ww
@@ -1062,11 +1064,13 @@ do while (item < nitems)
   end select
 end do
 
-print "(2a)", "Polarizabilities from file ", datafile
-unit=find_io(17)
-open(unit,file=datafile,status="old",iostat=ok)
-if (ok>0) call die ("Can't open file "//trim(datafile),.true.)
-call stream(unit)
+if (datafile .ne. "") then
+  print "(2a)", "Polarizabilities from file ", datafile
+  unit=find_io(17)
+  open(unit,file=datafile,status="old",iostat=ok)
+  if (ok>0) call die ("Can't open file "//trim(datafile),.true.)
+  call stream(unit)
+end if
 
 !  Read polarizability matrices individually
 do
@@ -2145,7 +2149,7 @@ INTEGER, INTENT(IN) :: rank
 !  rank=1 to include polarizabilities up to dipole
 !  rank=0 to include charge-flows only.
 
-DOUBLE PRECISION :: alphadd(3,3), w(1:4,1:4), d(1:4,1:4)
+DOUBLE PRECISION :: a(3,3), w(1:4,1:4), d(1:4,1:4), alpha_mean, gamma
 INTEGER :: i, j, ix, ka, kb
 
 !  Polarizability components are in multipole order (0,z,x,y) while
@@ -2153,7 +2157,7 @@ INTEGER :: i, j, ix, ka, kb
 INTEGER :: map(3)=[3,4,2]
 
 call axes(mol,.false.)
-alphadd = 0d0
+a = 0d0
 ka = firstp(mol)
 do while (ka > 0)
   kb = firstp(mol)
@@ -2167,9 +2171,9 @@ do while (ka > 0)
       w = matmul(w,transpose(d))
       do i = 1,3
         do j = 1,3
-          alphadd(i,j) = alphadd(i,j) + sx(i,ka)*sx(j,kb)*w(1,1)
+          a(i,j) = a(i,j) + sx(i,ka)*sx(j,kb)*w(1,1)
           if (rank > 0) then
-            alphadd(i,j) = alphadd(i,j)                                 &
+            a(i,j) = a(i,j)                                 &
                 + sx(i,ka)*w(1,map(j)) + w(map(i),1)*sx(j,kb)           &
                 + w(map(i),map(j))
           end if
@@ -2186,7 +2190,31 @@ print "(/2a)", "Total dipole-dipole polarizability for molecule ",     &
 if (rank == 0) print "(a)", "including charge-flow polarizabilities only"
 print "(a/(a,3f15.8))",                                                &
     "          x              y              z",                       &
-    " x ", alphadd(1,:), " y ", alphadd(2,:), " z ", alphadd(3,:)
+    " x ", a(1,:), " y ", a(2,:), " z ", a(3,:)
+
+alpha_mean = (a(1,1) + a(2,2) + a(3,3))/3d0
+
+if (a(1,2) == 0 .and. a(2,3) == 0 .and. a(3,1) == 0   &
+    .and. (abs(a(1,1) - a(2,2)) < 1d-6                          &
+    .or. abs(a(2,2) - a(3,3)) < 1d-6                            & 
+    .or. abs(a(3,3) - a(1,1)) < 1d-6)) then
+  !  Linear or symmetric top
+  if (abs(a(1,1) - a(2,2)) < 1d-6 .and. abs(a(2,2) - a(3,3)) < 1d-6) then
+    gamma = 0d0
+  else if (abs(a(1,1) - a(2,2)) < 1d-6) then
+    gamma = a(3,3) - (a(1,1) + a(2,2))/2d0
+  else if (abs(a(2,2) - a(3,3)) < 1d-6) then
+    gamma = a(1,1) - (a(2,2) + a(3,3))/2d0
+  else
+    gamma = a(2,2) - (a(3,3) + a(1,1))/2d0
+  end if
+else
+  gamma = sqrt(3d0*(a(1,2)**2 + a(2,3)**2 + a(3,1)**2)            &
+      + 0.5d0 * ((a(1,1)-a(2,2))**2 + (a(2,2)-a(3,3))**2 + (a(3,3)-a(1,1))**2))
+end if
+
+print "(a,f10.5,a,f10.5,a)", "Mean polarizability = ", alpha_mean,      &
+    "  anisotropy = ", gamma, " a.u."
 
 END SUBROUTINE total_dd