Fixed some bugs in force calculation (more to come)

The here fixed bugs are due to a missing division by 2*Pi wherever
the old cotra1 subroutine was replaced by a call to MATMUL.

force/force_a12.f90

@@ -10,6 +10,7 @@ CONTAINS
        we,jsp,ne,usdus,acof,bcof,e1cof,e2cof,&
        acoflo,bcoflo, results,f_a12)
     USE m_types
+    USE m_constants
     IMPLICIT NONE
 
     TYPE(t_results),INTENT(INOUT)   :: results
@@ -166,7 +167,7 @@ CONTAINS
 
              vec(:) = REAL(gv(:)) /atoms%neq(n)
 
-             gvint=MATMUL(cell%bmat,vec)
+             gvint=MATMUL(cell%bmat,vec)/tpi_const
              !
 
              vecsum(:) = zero

force/force_a21.F90

@@ -27,6 +27,7 @@ CONTAINS
     USE m_forcea21U
     USE m_tlmplm_store
     USE m_types
+    USE m_constants
     IMPLICIT NONE
 
     TYPE(t_results),INTENT(INOUT):: results
@@ -270,8 +271,8 @@ CONTAINS
              vec(:) = a21(:,natrun)
              vec2(:) = b4(:,natrun)
 
-             gvint=matmul(cell%bmat,vec)
-             gvint2=matmul(cell%bmat,vec2)
+             gvint=matmul(cell%bmat,vec)/tpi_const
+             gvint2=matmul(cell%bmat,vec2)/tpi_const
              vecsum(:) = zero
              vecsum2(:) = zero
 

force/geo.f90

@@ -39,6 +39,7 @@ CONTAINS
     USE m_bfgs
     USE m_bfgs0
     USE m_types
+    USE m_constants
     USE m_rinpXML
     USE m_winpXML
     IMPLICIT NONE
@@ -140,7 +141,7 @@ CONTAINS
     ELSE
        na = 0
        DO itype=1,atoms%ntype
-          tau0_i(:,itype)=MATMUL(cell%bmat,tau0(:,itype))
+          tau0_i(:,itype)=MATMUL(cell%bmat,tau0(:,itype))/tpi_const
           DO ieq = 1,atoms%neq(itype)
              na = na + 1
              jop = sym%invtab(atoms%ngopr(na))