! fortran_dialect=elf
!-----------------------------------------------------------------------
! PSCF - Polymer Self-Consistent Field Theory
! Copyright (2002-2016) Regents of the University of Minnesota
! contact: David Morse, morse012@umn.edu
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation. A copy of this license is included in
! the LICENSE file in the top-level PSCF directory. 
!-----------------------------------------------------------------------
!****m* scf/deform_mod
! MODULE
!   deform_mod 
! PURPOSE
!   Routines to treat prescribed deformation of unit cell. These may
!   be used calculate elastic properties in response to deformations
!   that lower the symmetry of the unperturbed crystal.
! AUTHOR
!   Chris Tyler (2002-2004)
! SOURCE
!-----------------------------------------------------------------------
module deform_mod 
   use const_mod  ! Defines constant long, and variable dim=1,2, or 3
   use io_mod     ! Standard interfaces for parameter io
   use group_mod  ! Defines types and operations for crystal groups
   use grid_mod   ! Utilities for fft grids
   use basis_mod  ! Utilities to make symmetrized basis functions
   implicit none
 
   PRIVATE

   ! Public procedures 
   PUBLIC :: make_deformed_basis
   PUBLIC :: deformation_subgroup
   PUBLIC :: deform_R_basis
   !***

contains

   subroutine make_deformed_basis( &
             R_basis,    &! real basis vectors, undeformed group
             G_basis,    &! reciprocal basis vectors, undeformed group
             group_name, &! file for undeformed space group
             deformation,&! deformation
             N_grids     &! maximum value of |G|
                      )
   !--------------------------------------------------------------------
   ! Analog of make_basis for group that is generated by finding the
   ! subgroup of the undeformed group that commutes with a deformation
   !--------------------------------------------------------------------
   use const_mod
   use group_mod
   use unit_cell_mod, only : make_G_basis
   use space_groups_mod
   implicit none

   real(long)  , intent(INOUT) :: R_basis(:,:)   ! real lattice basis 
   real(long)  , intent(OUT)   :: G_basis(:,:)   ! Reciprocal basis 
   character(len=*), intent(INOUT) :: group_name ! name of group
   integer     , intent(IN)    :: N_grids(3)     ! # of grid points
!   real(long),   intent(INOUT) :: Gabs_max       ! max. value of |G|
   real(long), intent(IN), dimension(:,:) :: deformation

   ! Local Variables
   integer :: i, j, k, l
   logical :: keep_cancel

   character(len = 3) :: Nchar

   type(group_type) :: group
   type(group_type) :: subgroup

   keep_cancel = .false.

   ! Make reciprocal lattice basis
   call make_G_basis(R_basis, G_basis )

   ! Read in and complete group
   call space_groups( group_name, group)
   call make_group( group, R_basis, G_basis)

   ! Find subgroup under deformation
   call deformation_subgroup( deformation, group, subgroup )

   ! Complete the deformed subgroup
   call deform_R_basis( R_basis, deformation )
   call make_G_basis( R_basis, G_basis )
   call make_group( subgroup, R_basis, G_basis )
   
   ! Make waves
   call make_waves(          & 
                G_basis,     &! (dim,dim), basis for reciprocal lattice
                R_basis,     &! (dim,dim), basis for Bravais lattice
                N_grids,     &! maximum allowed |G|
                subgroup,    &! space group 
                keep_cancel = .false. &
                        )

   ! group waves into stars related by space group symmetries
   call make_stars(subgroup)

   ! Make FFT grid of values of k^2 for reciprocal lattice vectors
   call make_ksq(G_basis)

   end subroutine make_deformed_basis
   !==================================================================


   !------------------------------------------------------------------
   !****p deform_mod/deformation_subgroup
   ! SUBROUTINE
   !    deformation subgroup
   ! PURPOSE
   !    Finds subgroup of group that commutes with a deformation tensor
   ! SOURCE
   !------------------------------------------------------------------
   subroutine deformation_subgroup( & 
                       deformation, & ! deformation tensor
                       group,       & ! original group
                       subgroup     & ! subgroup on output
                                 )
   use const_mod
   use group_mod
   implicit none

   real(long), dimension(:,:), intent(IN) :: deformation (dim,dim)
   type(group_type), intent(IN)           :: group
   type(group_type), intent(OUT)          :: subgroup
   !***
   
   ! Local arguments
   real(long), dimension(dim,dim) :: A, B
   integer :: i, k
   logical, dimension(dim,dim) :: compare
   
   real(long), parameter :: zero = 1.0d-8
   
   subgroup%order = 0
   do i = 1,group%order

      ! check if deformation*symmetry = symmetry *deformation
      A = matmul( deformation, group%s(i)%m )
      B = matmul( group%s(i)%m, deformation )
      compare = ( abs(A-B) < zero )
   
      ! If symmetry commutes with deformation, add to subgroup
      if ( all(compare) ) then
         subgroup%order = subgroup%order + 1
         k = subgroup%order
         subgroup%s(k)%m = group%s(i)%m
         subgroup%s(k)%v = group%s(i)%v
      endif

   enddo

   end subroutine deformation_subgroup
   !==================================================================
 
 
   !------------------------------------------------------------------
   !****p deform_mod/deform_R_basis
   ! SUBROUTINE
   !    make_dGsq(dGsq,dGG)
   ! PURPOSE
   !    Finds subgroup of group that commutes with a deformation tensor
   ! SOURCE
   !------------------------------------------------------------------
   subroutine deform_R_basis( R_basis, deformation)
   !------------------------------------------------------------------
   ! deform a R_basis according to the deformation 
   !------------------------------------------------------------------
   real(long), intent(INOUT) :: R_basis(:,:)
   real(long), intent(IN)    :: deformation(:,:)
 
   ! Local Variables
   real(long), dimension(dim,dim) :: def_basis
   integer :: i, j, k, l
 
   def_basis = 0.0_long
   
   do i = 1,dim
      def_basis(:,i) = matmul(deformation,R_basis(:,i))
   enddo
 
   R_basis = def_basis
 
   end subroutine deform_R_basis
   !==================================================================

end module deform_mod
