Merge pull request #77 from OrderN/hotfix-1.0.6-pre

Hot-fix for force and stress bugs
OrderN · Sep 11, 2020 · 087c2a6 · 087c2a6
2 parents f5075b9 + 25eb9f4
commit 087c2a6
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 19 deletions.
diff --git a/src/force_module.f90 b/src/force_module.f90
@@ -1248,7 +1248,7 @@ subroutine pulay_force(p_force, KE_force, fixed_potential, vary_mu,  &
                      i = bundle%ig_prim(iprim)
                      do isf = 1, natomf_species(bundle%species(iprim))
                         ! only accumulate phi-pulay force 3 times in total (not 9)
-                        if (flag_full_stress) then
+                        if (flag_stress .and. flag_full_stress) then
                            if (dir1 == dir2) then
                               p_force(dir1, i) = p_force(dir1, i) - &
                                  return_matrix_value(mat_tmp, np, ni, 0, 0, isf, isf, 1)
@@ -3150,7 +3150,6 @@ end subroutine matrix_diagonal
   subroutine matrix_diagonal_stress(matA, matB, stress, range, inode, diagonal)
 
     use datatypes
-    use global_module, only: atomic_stress
     use primary_module, only : bundle
     use matrix_module, only: matrix, matrix_halo
     use matrix_data, only : mat, halo

diff --git a/src/ion_electrostatic_module.f90 b/src/ion_electrostatic_module.f90
@@ -621,6 +621,10 @@ subroutine set_ewald(inode,ionode)
          ewald_g_vector_y(number_of_g_vectors), &
          ewald_g_vector_z(number_of_g_vectors),&
          ewald_g_factor(number_of_g_vectors),STAT=stat)
+    ewald_g_vector_x = zero
+    ewald_g_vector_y = zero
+    ewald_g_vector_z = zero
+    ewald_g_factor = zero
     if(stat/=0) &
          call cq_abort("set_ewald: error allocating ewald_g_vector ",&
          number_of_g_vectors,stat)
@@ -949,6 +953,8 @@ end subroutine structure_factor
 !!    Modifications to compute off-diagonal stress tensor elements
 !!   2019/05/08 zamaan
 !!    Added atomic stress contributions
+!!   2020/09/08 17:18 dave
+!!    Bug fix for reciprocal-space component of stress
 !!  SOURCE
 !!  
   subroutine ewald()
@@ -978,7 +984,7 @@ subroutine ewald()
     integer :: i, ip, ig_atom_beg, j, n, nc, ni, nj, nn, stat, direction, dir1, dir2, iatom, jatom
     real(double) :: arg_1, arg_2, coeff_1, dummy, ewald_real_energy, ewald_real_sum_inter, ewald_real_sum_intra, &
          &ewald_real_sum_ip, ewald_recip_energy, ewald_total_energy, g_dot_r, q_i, q_j, rij, rijx, rijy, rijz, &
-         rij_squared, x, y, z, one_over_g_squared
+         rij_squared, x, y, z, one_over_g_squared, delta_ab, one_over_four_gamma
     real(double), dimension(3) :: egv_n, rij_vec
     real(double), allocatable, dimension(:) :: ewald_recip_force_x, ewald_recip_force_y, ewald_recip_force_z, &
          &ewald_intra_force_x, ewald_intra_force_y, ewald_intra_force_z, &
@@ -1026,6 +1032,7 @@ subroutine ewald()
     call reg_alloc_mem(area_general,2*number_of_g_vectors,type_dbl)
     call structure_factor(inode,ionode)
     ! ... loop over reciprocal lattice vectors
+    one_over_four_gamma = one/(four*ewald_gamma)
     do n = 1, number_of_g_vectors
        ewald_recip_energy = ewald_recip_energy + &
             &(struc_fac_r(n)*struc_fac_r(n) + struc_fac_i(n)*struc_fac_i(n)) * &
@@ -1041,21 +1048,26 @@ subroutine ewald()
        one_over_g_squared = &
          one/(egv_n(1)*egv_n(1) + egv_n(2)*egv_n(2) + egv_n(3)*egv_n(3))
        if (flag_stress) then
-         do dir1=1,3
+          do dir1=1,3
            if (flag_full_stress) then
-             do dir2=1,3
-               ewald_recip_stress(dir1,dir2) = ewald_recip_stress(dir1,dir2)+ &
-                 ewald_g_factor(n) * (struc_fac_r(n)*struc_fac_r(n) + &
-                 struc_fac_i(n)*struc_fac_i(n)) * &
-                 ((two*egv_n(dir1)*egv_n(dir2) * & ! off-diagonal - zamaan
-                  (one/(four*ewald_gamma)) + one_over_g_squared - one))
+              do dir2=1,3
+                 if(dir2==dir1) then
+                    delta_ab = one
+                 else
+                    delta_ab = zero
+                 end if
+                 ewald_recip_stress(dir1,dir2) = ewald_recip_stress(dir1,dir2)+ &
+                      ewald_g_factor(n) * (struc_fac_r(n)*struc_fac_r(n) + &
+                      struc_fac_i(n)*struc_fac_i(n)) * &
+                      (two*egv_n(dir1)*egv_n(dir2) * &
+                      (one_over_four_gamma + one_over_g_squared) - delta_ab)
              end do
            else
              ewald_recip_stress(dir1,dir1) = ewald_recip_stress(dir1,dir1) + &
-               ewald_g_factor(n) * (struc_fac_r(n)*struc_fac_r(n) + &
-               struc_fac_i(n)*struc_fac_i(n)) * &
-               ((two*egv_n(dir1)*egv_n(dir1) * &
-                (one/(four*ewald_gamma)) + one_over_g_squared - one))
+                  ewald_g_factor(n) * (struc_fac_r(n)*struc_fac_r(n) + &
+                  struc_fac_i(n)*struc_fac_i(n)) * &
+                  (two*egv_n(dir1)*egv_n(dir1) * &
+                  (one_over_four_gamma + one_over_g_squared) - one)
            end if
          end do
        end if
@@ -1098,7 +1110,7 @@ subroutine ewald()
     enddo
     ! Correctly scale stress
     if (flag_stress) then
-      do dir1=1,3 
+       do dir1=1,3 
         do dir2=1,3
           ewald_recip_stress(dir1,dir2) = ewald_recip_stress(dir1,dir2) * &
                                           two * pi/(ewald_real_cell_volume)
@@ -1362,7 +1374,7 @@ subroutine ewald()
     call gsum(ewald_total_force_z,ni_in_cell)
     ! SYM 2014/10/22 14:34: Summ all the stress contributions
     if (flag_stress) then
-      do dir1=1,3
+       do dir1=1,3
          do dir2=1,3
             ion_interaction_stress(dir1,dir2) = &
               ewald_intra_stress(dir1,dir2) + ewald_inter_stress(dir1,dir2) + &

diff --git a/src/pseudo_tm_module.f90 b/src/pseudo_tm_module.f90
@@ -957,6 +957,8 @@ end subroutine set_tm_pseudo
 !!    Off-diagonal stress tensor elements
 !!   2019/05/08 zamaan
 !!    Added atomic stress contributions
+!!   2020/09/08 17:18 dave (via tsuyoshi)
+!!    Bug fix for local_HF_force
 !!  SOURCE
 !!
   subroutine loc_pp_derivative_tm ( hf_force, density, size )
@@ -1283,7 +1285,7 @@ subroutine loc_pp_derivative_tm ( hf_force, density, size )
                               fr_2 = front * r(dir1)
                               HF_force(dir1,ig_atom) = &
                                 HF_force(dir1,ig_atom) + &
-                                fr_1(dir1) * elec_here * fr_2(dir1)
+                                fr_1(dir1) * elec_here + fr_2(dir1)
                               if (flag_stress) then
                                 if (flag_full_stress) then
                                   do dir2=1,3
@@ -1293,8 +1295,8 @@ subroutine loc_pp_derivative_tm ( hf_force, density, size )
                                       r(dir2) * r_from_i
                                     if (flag_atomic_stress) then
                                       atomic_stress(dir1,dir2,ig_atom) = &
-                                        atomic_stress(dir1,dir2,ig_atom) * &
-                                        (fr_1(dir1) * elec_here + fr_2(dir2))*&
+                                        atomic_stress(dir1,dir2,ig_atom) + &
+                                        (fr_1(dir1) * elec_here + fr_2(dir1))*&
                                         r(dir2) * r_from_i
                                     end if
                                   end do