Adding rotation calculation

src/noft/noft_module.F

@@ -54,6 +54,7 @@ module noft
 #include "noft_load_vectors.fi"
 #include "noft_local.fi"
 #include "noft_monte_carlo_range.fi"
+#include "noft_monte_carlo_rotation.fi"
 #include "noft_monte_carlo_step.fi"
 #include "noft_ortho_operators.fi"
 #include "noft_orthogonalize_vecs.fi"

src/noft/noft_monte_carlo_range.fi

@@ -158,13 +158,18 @@
       real(kind=dp) :: scale_b
       real(kind=dp) :: scale_c
       real(kind=dp) :: thresh
+      real(kind=dp) :: rangemx
 !
       scale_a = noft_parameters%monte_carlo%scale_a
       scale_b = noft_parameters%monte_carlo%scale_b
       scale_c = noft_parameters%monte_carlo%scale_c
       thresh  = noft_parameters%monte_carlo%step_threshold
+      rangemx = noft_parameters%monte_carlo%step_maximum
       call noft_access(ga_range,range,ilo,ihi)
       allocate(step(ilo(1):ihi(1),ilo(2):ihi(2)))
+      do ii = 1, ndim
+        ild(ii) = ihi(ii) - ilo(ii) + 1
+      enddo
       call nga_get(ga_step,ilo,ihi,step,ild)
       do jj = ilo(2), ihi(2)
         do ii = ilo(1), ihi(1)
@@ -174,6 +179,7 @@
      &                   + (scale_b*(2.0d0*step(ii,jj)-range(ii,jj))
      &                      + scale_c*range(ii,jj))*(1.0d0-scale_a)
             range(ii,jj) = max(range(ii,jj),0.1d0*thresh)
+            range(ii,jj) = min(range(ii,jj),rangemx)
           endif
         enddo
       enddo

src/noft/noft_monte_carlo_rotation.fi

@@ -0,0 +1,167 @@
+!-----------------------------------------------------------------------
+!
+!> \brief Create Monte Carlo rotation matrices
+!>
+!> This subroutine creates all the global arrays for a rotation data
+!> type.
+!>
+      subroutine noft_monte_carlo_rotation_create(noft_parameters,
+     &                                            noft_rotation)
+      implicit none
+#include "global.fh"
+#include "mafdecls.fh"
+#include "errquit.fh"
+      !> The calculation parameters
+      type(noft_parameter_tp),            intent(in)  :: noft_parameters
+      !> The rotation
+      type(noft_monte_carlo_rotation_tp), intent(out) :: noft_rotation
+!
+!     Local
+!
+      character(len=30), parameter :: pname
+     &   = "noft_monte_carlo_rotation_create: "
+      integer :: nmo
+      integer :: nto
+      integer :: mo_a, mo_b
+      integer :: to_a, to_b
+!
+      nmo = noft_parameters%nmo
+      nto = noft_parameters%nto
+!
+      if (.not.ga_create(MT_DBL,nmo,nmo,"rotation mo a",-1,-1,mo_a))
+     &  call errquit(pname//"failed to create MO_a",10,GA_ERR)
+      if (.not.ga_create(MT_DBL,nmo,nmo,"rotation mo b",-1,-1,mo_b))
+     &  call errquit(pname//"failed to create MO_b",20,GA_ERR)
+      if (.not.ga_create(MT_DBL,nto,nto,"rotation to a",-1,-1,to_a))
+     &  call errquit(pname//"failed to create TO_a",30,GA_ERR)
+      if (.not.ga_create(MT_DBL,nto,nto,"rotation mo b",-1,-1,to_b))
+     &  call errquit(pname//"failed to create TO_b",40,GA_ERR)
+!
+      noft_rotation%mo_a = mo_a
+      noft_rotation%mo_b = mo_b
+      noft_rotation%to_a = to_a
+      noft_rotation%to_b = to_b
+!
+      end subroutine noft_monte_carlo_rotation_create
+!
+!-----------------------------------------------------------------------
+!
+!> \brief Destroy Monte Carlo rotation matrices
+!>
+!> Simply clean up all the global arrays associated with a Monte Carlo
+!> rotation data type.
+!>
+      subroutine noft_monte_carlo_rotation_destroy(noft_rotation)
+      implicit none
+#include "global.fh"
+#include "errquit.fh"
+      !> The Monte Carlo rotation
+      type(noft_monte_carlo_rotation_tp), intent(inout) :: noft_rotation
+      character(len=30), parameter :: pname
+     &   = "noft_monte_carlo_rotation_destroy: "
+!
+      if (.not.ga_destroy(noft_rotation%mo_a))
+     &  call errquit(pname//"failed to destroy mo_a",10,GA_ERR)
+      if (.not.ga_destroy(noft_rotation%mo_b))
+     &  call errquit(pname//"failed to destroy mo_b",20,GA_ERR)
+      if (.not.ga_destroy(noft_rotation%to_a))
+     &  call errquit(pname//"failed to destroy to_a",30,GA_ERR)
+      if (.not.ga_destroy(noft_rotation%to_b))
+     &  call errquit(pname//"failed to destroy to_b",40,GA_ERR)
+!
+      end subroutine noft_monte_carlo_rotation_destroy
+!
+!-----------------------------------------------------------------------
+!
+!> \brief Compute Monte Carlo rotation matrix
+!>
+!> Given a skew symmetric step matrix in terms of angles compute
+!> the corresponding rotation matrix.
+!>
+      subroutine noft_monte_carlo_compute_rotation_1(ga_rotation,
+     &                                               ga_step)
+      implicit none
+#include "errquit.fh"
+#include "global.fh"
+      !> The rotation matrix
+      integer, intent(inout) :: ga_rotation
+      !> The skew symmetric step
+      integer, intent(in)    :: ga_step
+!
+!     Local
+!
+      character(len=30), parameter :: pname
+     &  = "noft_monte_carlo_compute_rotation_1: "
+      integer, parameter :: mdim = 2
+      integer :: idim(mdim)
+      integer :: ilo(mdim)
+      integer :: ihi(mdim)
+      integer :: ild(mdim)
+      integer :: ndim
+      integer :: iproc
+      integer :: nproc
+      integer :: itype
+      integer :: ncol ! the number of columns for this processor
+      integer :: i1, i2
+      !real(kind=dp) :: cc
+      real(kind=dp) :: cd
+      !real(kind=dp) :: cs
+      real(kind=dp), allocatable :: buf(:,:)
+!
+      nproc = ga_nnodes()
+      iproc = ga_nodeid()
+      call nga_inquire(ga_step,itype,ndim,idim)
+      if (ndim.eq.2) then
+        ilo(1) = 1
+        ihi(1) = idim(1)
+        ncol = (idim(2)+nproc-1)/nproc
+        ilo(2) = iproc*ncol+1
+        ihi(2) = min((iproc+1)*ncol,idim(2))
+        ild(1) = ihi(1)-ilo(1)+1
+        ild(2) = ihi(2)-ilo(2)+1
+        if (ilo(2).le.ihi(2)) then
+          allocate(buf(ilo(1):ihi(1),ilo(2):ihi(2)))
+          call nga_get(ga_step,ilo,ihi,buf,ild)
+          do i2 = ilo(2), ihi(2)
+            cd = 1.0_dp
+            do i1 = ilo(1), ihi(1)
+              cd = cd * cos(buf(i1,i2))
+              buf(i1,i2) = sin(buf(i1,i2))
+            enddo
+            buf(i2,i2) = cd
+          enddo
+          call nga_put(ga_rotation,ilo,ihi,buf,ild)
+          deallocate(buf)
+        endif
+      else
+        call errquit(pname//"invalid number of dimensions",
+     &               10, CAPMIS_ERR)
+      endif
+!
+      end subroutine noft_monte_carlo_compute_rotation_1
+!
+!-----------------------------------------------------------------------
+!
+      subroutine noft_monte_carlo_compute_rotation(noft_rotation,
+     &                                             noft_step)
+      implicit none
+      !> The rotation
+      type(noft_monte_carlo_rotation_tp), intent(inout) :: noft_rotation
+      !> The step
+      type(noft_monte_carlo_step_tp),     intent(in)    :: noft_step
+
+!
+      call ga_sync()
+      call noft_monte_carlo_compute_rotation_1(
+     &         noft_rotation%mo_a,noft_step%mo_a)
+      call noft_monte_carlo_compute_rotation_1(
+     &         noft_rotation%mo_b,noft_step%mo_b)
+      call noft_monte_carlo_compute_rotation_1(
+     &         noft_rotation%to_a,noft_step%to_a)
+      call noft_monte_carlo_compute_rotation_1(
+     &         noft_rotation%to_b,noft_step%to_b)
+      call ga_sync()
+!
+      end subroutine noft_monte_carlo_compute_rotation
+!
+!-----------------------------------------------------------------------

src/noft/noft_monte_carlo_step.fi

@@ -73,6 +73,24 @@
 !
 !-----------------------------------------------------------------------
 !
+!> \brief Choose a random step for one set of coordinates
+!>
+!> The step is given in terms of angles. The rotation that we will
+!> construct from this is given by
+!>
+!>      (R_ii R_ij)   ( c  -s )
+!>  R = (         ) = (       )
+!>      (R_ji R_jj)   ( s   c )
+!>
+!> where R_ij = -R_ji (see e.g.
+!> https://en.wikipedia.org/wiki/Jacobi_rotation).
+!> If the step is given in terms of angles then \f$c = \cos(\theta)\f$
+!> and \f$s = \sin(\theta)\f$. Note that \f$\cos\f$ is symmetric around
+!> \f$0\f$ whereas \f$\sin\f$ is anti-symmetric. That means that if 
+!> we ensure that the step matrix is skew symmetric then the skew
+!> symmetry of the rotation matrix is automatically obtained.
+!> Therefore we choose the step matrix so that it is skew symmtric.
+!>
       subroutine noft_monte_carlo_step_choose_1(ga_step,ga_range)
       implicit none
       !> The Monte Carlo step global array handle
@@ -88,9 +106,14 @@
       integer :: ild(ndim)
       integer :: ii, jj
       real(kind=dp), pointer     :: step(:,:)
+      real(kind=dp), allocatable :: stept(:,:)
       real(kind=dp), allocatable :: range(:,:)
+      real(kind=dp) :: onem
 !
+      onem = -1.0d0
+      call ga_zero(ga_step)
       call noft_access(ga_step,step,ilo,ihi)
+      allocate(stept(ilo(2):ihi(2),ilo(1):ihi(1)))
       allocate(range(ilo(1):ihi(1),ilo(2):ihi(2)))
       do ii = 1, ndim
         ild(ii) = ihi(ii) - ilo(ii) + 1
@@ -99,11 +122,15 @@
       call random_number(step)
       do jj = ilo(2), ihi(2)
         do ii = ilo(1), ihi(1)
-          step(ii,jj) = 2.0d0*step(ii,jj)*range(ii,jj) - range(ii,jj)
+          step(ii,jj)  = 2.0d0*step(ii,jj)*range(ii,jj) - range(ii,jj)
+          stept(jj,ii) = step(ii,jj)
         enddo
       enddo
       deallocate(range)
       call noft_release(ga_step,step,ilo,ihi)
+      call ga_acc(ga_step,ilo(2),ihi(2),ilo(1),ihi(1),stept,ild(2),onem)
+      call ga_sync()
+      deallocate(stept)
 !
       end subroutine noft_monte_carlo_step_choose_1
 !