Long-range should be be properly symmetric now

src/extract_forceconstants/main.f90

@@ -51,7 +51,6 @@ program extract_forceconstants
 ! We need a forcemap to continue.
 getfrcmap: block
     type(lo_interaction_tensors) :: slt
-    real(r8) :: t0
 
     ! Either read it from file or create it
     if (opts%readforcemap) then
@@ -62,8 +61,8 @@ program extract_forceconstants
         end if
         call map%read_from_hdf5(uc, 'infile.forcemap.hdf5', opts%verbosity)
         ! Update constraints that depend on structure
-        call lo_secondorder_rot_herm_huang(map, uc, map%constraints%eq2, map%constraints%neq2, &
-                                           opts%rotationalconstraints, opts%huanginvariance, opts%hermitian)
+        !call lo_secondorder_rot_herm_huang(map, uc, map%constraints%eq2, map%constraints%neq2, &
+        !                                   opts%rotationalconstraints, opts%huanginvariance, opts%hermitian)
     else
         ! Now we have to calculate the whole thing.
         call uc%classify('wedge', timereversal=.true.)
@@ -90,12 +89,6 @@ program extract_forceconstants
         call map%generate(uc, ss, polarcorrectiontype=opts%polarcorrectiontype, &
                           st=slt, mw=mw, mem=mem, verbosity=opts%verbosity, devmode=opts%devmode)
 
-        ! Create the constraints
-        t0 = walltime()
-        call map%forceconstant_constraints(uc, opts%rotationalconstraints, &
-                                           opts%huanginvariance, opts%hermitian, opts%verbosity + 10)
-        if (mw%talk) write (*, *) '... got constraints in ', tochar(walltime() - t0), 's'
-        ! Maybe print forcemap to file.
         if (opts%printforcemap .and. mw%talk) then
             call map%write_to_hdf5('outfile.forcemap.hdf5', opts%verbosity)
         end if
@@ -339,14 +332,22 @@ program extract_forceconstants
     end if
     if (map%have_fc_pair) then
         call map%get_secondorder_forceconstant(uc, fc2, mem, opts%verbosity)
+        call fc2%get_elastic_constants(uc, mw, opts%verbosity)
         if (mw%talk) then
             call fc2%writetofile(uc, 'outfile.forceconstant')
-            call fc2%get_elastic_constants(uc)
             write (*, *) ''
-            write (*, *) 'ELASTIC CONSTANTS (GPa):'
+            write (*, *) 'ELASTIC CONSTANTS (short range) (GPa):'
             do i = 1, 6
                 write (*, "(6(3X,F15.5))") lo_chop(fc2%elastic_constants_voigt(:, i)*lo_pressure_HartreeBohr_to_GPa, lo_tol)
             end do
+            write (*, *) 'ELASTIC CONSTANTS (long range) (GPa):'
+            do i = 1, 6
+                write (*, "(6(3X,F15.5))") lo_chop(fc2%elastic_constants_voigt_longrange(:, i)*lo_pressure_HartreeBohr_to_GPa, lo_tol)
+            end do
+            write (*, *) 'ELASTIC CONSTANTS (total) (GPa):'
+            do i = 1, 6
+                write (*, "(6(3X,F15.5))") lo_chop((fc2%elastic_constants_voigt(:, i) + fc2%elastic_constants_voigt_longrange(:, i))*lo_pressure_HartreeBohr_to_GPa, lo_tol)
+            end do
         end if
     end if
     if (map%have_fc_triplet) then

src/extract_forceconstants/options.f90

@@ -124,17 +124,17 @@ subroutine parse(opts)
                       &These can be used to find the finite temperature equilibrium structure.', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--readforcemap',hidden=.true., &
+    call cli%add(switch='--readforcemap', hidden=.true., &
                  help='Read `infile.forcemap.hdf5` from file instead of calculating &
                       &all symmetry relations. Useful for sets of calculations with the same structure.', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--readirreducible',hidden=.true., &
+    call cli%add(switch='--readirreducible', hidden=.true., &
                  help='Read the irreducible forceconstants from `infile.irrifc_*` &
                       &instead of solving for them. This option requires an `infile.forcemap.hdf5`, as above.', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--potential_energy_differences', switch_ab='-U0',hidden=.true., &
+    call cli%add(switch='--potential_energy_differences', switch_ab='-U0', hidden=.true., &
                  help='Calculate the difference in potential energy from the &
                       &simulation and the forceconstants to determine U0.', &
                  help_markdown='As referenced in the thermodynamics section of &
@@ -146,12 +146,12 @@ subroutine parse(opts)
                                &This number should be added to the appropriate phonon free energy.', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--printforcemap',hidden=.true., &
+    call cli%add(switch='--printforcemap', hidden=.true., &
                  help='Print `outfile.forcemap.hdf5` for reuse.', &
                  required=.false., act='store_true', def='.false.', error=lo_status)
     if (lo_status .ne. 0) stop
-    call cli%add(switch='--temperature',&
-                 help='Temperature for self-consistent solver.',&
+    call cli%add(switch='--temperature', &
+                 help='Temperature for self-consistent solver.', &
                  required=.false., act='store', def='-1', error=lo_status)
     if (lo_status .ne. 0) stop
 

src/libolle/Makefile

@@ -51,7 +51,6 @@ $(OBJECT_PATH)lo_symmetry_of_interactions_helpers.o \
 $(OBJECT_PATH)lo_symmetry_of_interactions_tuplets.o \
 $(OBJECT_PATH)lo_symmetry_of_interactions_nullspace.o \
 $(OBJECT_PATH)lo_longrange_electrostatics.o \
-$(OBJECT_PATH)lo_longrange_electrostatics_optz.o \
 $(OBJECT_PATH)lo_longrange_electrostatics_dynmat.o \
 $(OBJECT_PATH)lo_dielectric_interaction.o \
 $(OBJECT_PATH)lo_dielectric_interaction_matrixelements.o \
@@ -299,11 +298,6 @@ $(OBJECT_PATH)lo_longrange_electrostatics.o:\
 	$(OBJECT_PATH)type_crystalstructure.o \
     $(OBJECT_PATH)gottochblandat.o
 	$(FC) $(FFLAGS) -c lo_longrange_electrostatics.f90 -o $@
-$(OBJECT_PATH)lo_longrange_electrostatics_optz.o:\
-	lo_longrange_electrostatics_optz.f90 \
-    $(OBJECT_PATH)lo_longrange_electrostatics.o \
-	$(OBJECT_PATH)quadratures_stencils.o
-	$(FC) $(FFLAGS) -c lo_longrange_electrostatics_optz.f90 -o $@
 $(OBJECT_PATH)lo_longrange_electrostatics_dynmat.o:\
 	lo_longrange_electrostatics_dynmat.f90 \
     $(OBJECT_PATH)lo_longrange_electrostatics.o

src/libolle/gottochblandat.f90

@@ -468,11 +468,12 @@ module pure subroutine lo_transpositionmatrix(tm)
 #endif
         real(flyt), dimension(:,:), intent(out) :: tm
     end subroutine
-    module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tolerance,subset,weights,gramified)
+    module subroutine lo_linear_least_squares(A,B,x,constraint_C,constraint_D,nconstraints,tolerance,subset,weights,gramified)
         real(flyt), dimension(:,:), intent(in) :: A
         real(flyt), dimension(:), intent(in) :: B
         real(flyt), dimension(:), intent(out) :: x
-        real(flyt), dimension(:,:), intent(in), optional :: zeroconstraints
+        real(flyt), dimension(:,:), intent(in), optional :: constraint_C
+        real(flyt), dimension(:), intent(in), optional :: constraint_D
         integer, intent(in), optional :: nconstraints
         real(flyt), intent(in), optional :: tolerance
         integer, intent(in), dimension(:), optional :: subset
@@ -1540,10 +1541,6 @@ subroutine lo_stop_gracefully(msg,exitcode,filename,line)
         write(*,*) 'exit code 6: I/O error'
     case(7)
         write(*,*) 'exit code 7: MPI error'
-    case(8)
-        write(*,*) 'exit code 8: Feature removed'
-    case(9)
-        write(*,*) 'exit code 9: Bad inputs.'
     end select
     write(*,*) ''
     do i=1,size(msg)

src/libolle/gottochblandat_linalg.f90

@@ -1411,15 +1411,17 @@ pure module subroutine lo_transpositionmatrix(tm)
 end subroutine
 
 !> wrapper around LAPACK so that it is not destructive. Solves Ax=B, perhaps with constraints such that Cx=0, with a bunch of fancy options.
-module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tolerance,subset,weights,gramified)
+module subroutine lo_linear_least_squares(A,B,x,constraint_C,constraint_D,nconstraints,tolerance,subset,weights,gramified)
     !> matrix A
     real(flyt), dimension(:,:), intent(in) :: A
     !> vector B
     real(flyt), dimension(:), intent(in) :: B
     !> solution
     real(flyt), dimension(:), intent(out) :: x
     !> linear constraints
-    real(flyt), dimension(:,:), intent(in), optional :: zeroconstraints
+    real(flyt), dimension(:,:), intent(in), optional :: constraint_C
+    !> D in Cx=D
+    real(flyt), dimension(:), intent(in), optional :: constraint_D
     !> number of constraints
     integer, intent(in), optional :: nconstraints
     !> specify a tolerance
@@ -1447,11 +1449,11 @@ module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tol
         ! least more robust than just a plain application.
 
         ! Initial sanity test for the constraints
-        if ( present(zeroconstraints) ) then
+        if ( present(constraint_C) ) then
             if ( .not.present(nconstraints) ) then
                 call lo_stop_gracefully(['Need to provide the number of constraints to the llsq'],lo_exitcode_param,__FILE__,__LINE__)
             elseif ( nconstraints .gt. 0 ) then
-                if ( nconstraints .ne. size(zeroconstraints,1) ) then
+                if ( nconstraints .ne. size(constraint_C,1) ) then
                     call lo_stop_gracefully(['Inconsistent number of constraints'],lo_exitcode_param,__FILE__,__LINE__)
                 endif
             endif
@@ -1534,7 +1536,7 @@ module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tol
                 lo_allocate(vC(nc,size(xind)))
                 vC=0.0_flyt
                 do i=1,size(xind)
-                    vC(:,i)=zeroconstraints(:,xind(i))
+                    vC(:,i)=constraint_C(:,xind(i))
                 enddo
                 call lo_compress_equations(vC,l,wC,trans=.false.,tolerance=lo_sqtol)
                 if ( l .ne. nc ) then
@@ -1556,14 +1558,17 @@ module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tol
             ! first copy the response vector
             vB=0.0_flyt
             vB(1:nx)=B
+            do i=1,nc
+                vB(nx+i)=constraint_D(i)
+            enddo
             ! then the full thingy
             wA=0.0_flyt
             wA(1:nx,1:nx)=A
             do j=1,nx
             do i=1,nc
                 l=nx+i
-                wA(l,j)=zeroconstraints(i,j)
-                wA(j,l)=zeroconstraints(i,j)
+                wA(l,j)=constraint_C(i,j)
+                wA(j,l)=constraint_C(i,j)
             enddo
             enddo
         else
@@ -1572,7 +1577,7 @@ module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tol
             wA=A
             if ( nc .gt. 0 ) then
                 lo_allocate(wC(nc,size(A,2)))
-                wC=zeroconstraints
+                wC=constraint_C
             endif
             ! And a copy of the prediction vector
             lo_allocate(vB(size(B)))
@@ -1698,22 +1703,24 @@ module subroutine lo_linear_least_squares(A,B,x,zeroconstraints,nconstraints,tol
 
     ! Now for the final touch.
     finalize: block
-        real(flyt), dimension(:), allocatable :: viol
-
-        ! Check if the constraints are satisfied
-        if ( present(zeroconstraints) ) then
-        if ( nconstraints .gt. 0 ) then
-            lo_allocate(viol(size(zeroconstraints,2)))
-            viol=0.0_flyt
-            call dgemv('N',size(zeroconstraints,1),size(zeroconstraints,2),1.0_flyt,&
-                       zeroconstraints,size(zeroconstraints,1),X,1,0.0_flyt,viol,1)
-            ! Check if they are not satisfied. If so, fix it.
-            if ( sum(abs(viol)) .gt. tol ) then
-                call lo_enforce_linear_constraints(zeroconstraints,x,tol,nosquare=.true.)
-            endif
-            lo_deallocate(viol)
-        endif
-        endif
+        !real(flyt), dimension(:), allocatable :: viol
+
+        ! ! Check if the constraints are satisfied
+        ! if ( present(zeroconstraints) ) then
+        ! if ( nconstraints .gt. 0 ) then
+        !     lo_allocate(viol(size(zeroconstraints,2)))
+        !     viol=0.0_flyt
+        !     ! call dgemv('N',size(zeroconstraints,1),size(zeroconstraints,2),1.0_flyt,&
+        !     !            zeroconstraints,size(zeroconstraints,1),X,1,0.0_flyt,viol,1)
+        !     ! Check if they are not satisfied. If so, fix it.
+        !     if ( sum(abs(viol)) .gt. tol ) then
+        !         write(*,*) 'Constraints not satisfied!'
+        !         stop
+        !         !call lo_enforce_linear_constraints(zeroconstraints,x,tol,nosquare=.true.)
+        !     endif
+        !     lo_deallocate(viol)
+        ! endif
+        ! endif
         ! Chop away small numbers
         X=lo_chop(X,tol*1E-2_flyt)
     end block finalize

src/libolle/ifc_solvers.f90

@@ -15,7 +15,7 @@ module ifc_solvers
 use lo_memtracker, only: lo_mem_helper
 use type_crystalstructure, only: lo_crystalstructure
 use type_qpointmesh, only: lo_qpoint
-use type_forcemap, only: lo_forcemap
+use type_forcemap, only: lo_forcemap,lo_coeffmatrix_unitcell_Z_singlet
 use type_mdsim, only: lo_mdsim
 use type_blas_lapack_wrappers, only: lo_gemv
 use lo_longrange_electrostatics, only: lo_ewald_parameters
@@ -304,14 +304,15 @@ module subroutine report_diagnostics(map, tp, ih, dh, mw, mem, verbosity)
 end interface
 ! ifc_solvers_fastugly
 interface
-    module subroutine lo_solve_for_irreducible_ifc_fastugly(map, uc, ss, sim, mw, mem, verbosity)
+    module subroutine lo_solve_for_irreducible_ifc_fastugly(map, uc, ss, sim, mw, mem, verbosity, fix_secondorder)
         type(lo_forcemap), intent(inout) :: map
         type(lo_crystalstructure), intent(in) :: uc
         type(lo_crystalstructure), intent(in) :: ss
         type(lo_mdsim), intent(in) :: sim
         type(lo_mpi_helper), intent(inout) :: mw
         type(lo_mem_helper), intent(inout) :: mem
         integer, intent(in) :: verbosity
+        logical, intent(in) :: fix_secondorder
     end subroutine
 end interface
 
@@ -462,6 +463,107 @@ subroutine lo_solve_for_irreducible_ifc(map, sim, uc, ss, solver, mw, mem, verbo
     end block dielectric
     call mem%tock(__FILE__, __LINE__, mw%comm)
 
+    ! As soon as we know the Born charges we can construct the IFC constraints
+    buildconstraints: block
+        type(lo_ewald_parameters) :: ew
+        complex(r8), dimension(:,:,:,:), allocatable :: D0
+        real(r8), dimension(:,:,:,:), allocatable :: rottensor
+        real(r8), dimension(:,:,:), allocatable :: wZ,wD
+        real(r8), dimension(:,:), allocatable :: coeff_Z
+        real(r8), dimension(:), allocatable :: v0
+        real(r8), dimension(:), allocatable :: hermitian_rhs, huang_rhs, rotational_rhs
+        real(r8), dimension(3,3,3,3) :: bracket
+        real(r8), dimension(3,3) :: eps,m0,m1
+        integer :: i,j,k,l,ii
+        integer :: a1,a2
+
+        allocate(hermitian_rhs(9*uc%na))
+        allocate(rotational_rhs(27*uc%na))
+        allocate(huang_rhs(81))
+        hermitian_rhs=0.0_r8
+        rotational_rhs=0.0_r8
+        huang_rhs=0.0_r8
+
+        if ( map%polar .gt. 0 ) then
+            ! First thing we need is the polar dynamical matrix to know
+            ! how we should build the Hermiticity constraints.
+            eps=lo_unflatten_2tensor(matmul(map%eps_global_shell%coeff, map%xuc%x_eps_global))
+
+            allocate(wZ(3,3,uc%na))
+            allocate(wD(3,3,uc%na))
+            allocate(v0(9*uc%na))
+            allocate(coeff_Z(uc%na*9,map%xuc%nx_Z_singlet))
+            allocate(D0(3,3,uc%na,uc%na))
+            wZ=0.0_r8
+            wD=0.0_r8
+            v0=0.0_r8
+            coeff_Z=0.0_r8
+            D0=0.0_r8
+
+            call lo_coeffmatrix_unitcell_Z_singlet(map, coeff_Z)
+            v0=matmul(coeff_Z,map%xuc%x_Z_singlet)
+            do a1 = 1, uc%na
+                wZ(:, :, a1) = lo_unflatten_2tensor(v0((a1 - 1)*9 + 1:a1*9))
+            end do
+
+            call ew%set(uc, eps, 2, 1E-20_r8, verbosity=-1)
+            call ew%longrange_dynamical_matrix(uc, [0.0_r8, 0.0_r8, 0.0_r8], wZ, wD, eps, D0, reconly=.true., chgmult=.true.)
+
+            do a1=1,uc%na
+                m0=0.0_r8
+                m1=0.0_r8
+                do a2=1,uc%na
+                    m0=m0+real(d0(:,:,a1,a2),r8)
+                    m1=m1+real(d0(:,:,a2,a1),r8)
+                enddo
+                m0=m0-m1
+                hermitian_rhs( (a1-1)*9+1:a1*9 ) = -lo_flattentensor(m0)
+            enddo
+            hermitian_rhs = lo_chop(hermitian_rhs,1E-12_r8)
+
+            ! Next up we do the huang + rotational invariances
+            allocate(rottensor(3,3,3,uc%na))
+            rottensor=0.0_r8
+            call ew%longrange_elastic_constant_bracket( uc, wZ, eps, bracket, rottensor, reconly=.true., mw=mw, verbosity=verbosity)
+
+            huang_rhs=0.0_r8
+            ii=0
+            do i=1,3
+            do j=1,3
+            do k=1,3
+            do l=1,3
+                ii=ii+1
+                huang_rhs(ii) = bracket(k,l,i,j) - bracket(i,j,k,l)
+            enddo
+            enddo
+            enddo
+            enddo
+            huang_rhs = lo_chop(huang_rhs,1E-12_r8)
+
+            ii=0
+            do a1=1,uc%na
+            do i=1,3
+            do j=1,3
+            do k=1,3
+                ii=ii+1
+                rotational_rhs(ii) = rottensor(i,j,k,a1)-rottensor(i,k,j,a1)
+            enddo
+            enddo
+            enddo
+            enddo
+            rotational_rhs = lo_chop(rotational_rhs,1E-12_r8)
+
+        else
+            ! We just have zeros
+            huang_rhs=0.0_r8
+            hermitian_rhs=0.0_r8
+            rotational_rhs=0.0_r8
+        endif
+
+        call map%forceconstant_constraints(uc,rotational=.true.,huanginvariances=.true.,hermitian=.true.,hermitian_rhs=hermitian_rhs,huang_rhs=huang_rhs,rotational_rhs=rotational_rhs,verbosity=verbosity)
+        !call map%forceconstant_constraints(uc,rotational=.false.,huanginvariances=.true.,hermitian=.true.,hermitian_rhs=hermitian_rhs,huang_rhs=huang_rhs,verbosity=verbosity)
+    end block buildconstraints
+
     ! Get some things that are always needed
     !call mem%tick()
     ifc: block