bug fix, freezing with gcc ver>8, resolves issue #12

Author: saroun

src/lattice.f

@@ -76,9 +76,8 @@ end subroutine LAT_CLEAR
       subroutine LAT_INIT(LAT)
 ! calculate dependent fields in LAT (metric tensor, cell volume, ...)
 !-------------------------------------------------------------
-      TYPE (TCR_LATTICE)  :: LAT
-      LAT%VOL=GET_CELLVOL(LAT,298.D0)
-      call GET_METRIC(LAT,LAT%G)
+      TYPE (TCR_LATTICE) :: LAT
+      call GET_METRIC(LAT)
       end subroutine LAT_INIT
 
 !--------------------------------------------------------
@@ -226,32 +225,34 @@ real(kind(1.D0)) function GET_CELLVOL(LAT,T)
       end function GET_CELLVOL
 
 !-------------------------------------------------------------
-      subroutine GET_METRIC(LAT,G)
+      subroutine GET_METRIC(LAT)
 ! calculate metric tensor for reciprocal lattice at T=298 K
 !-------------------------------------------------------------
-      TYPE (TCR_LATTICE),intent(in)  :: LAT
-      real(kind(1.D0)),intent(out) :: G(3,3)
+      TYPE (TCR_LATTICE) :: LAT
+      real(kind(1.D0)) :: G(3,3)
       integer :: k
-      real(kind(1.D0)) :: AX(3),SINA(3),COSA(3),VOL
-      VOL=GET_CELLVOL(LAT,298.D0)
-      if (VOL.gt.0.D0) then
+      real(kind(1.D0)) :: AX(3),SINA(3),COSA(3),VOL2
+	  LAT%VOL=GET_CELLVOL(LAT,298.D0)
+      if (LAT%VOL.gt.0.D0) then
+	    VOL2 = LAT%VOL**2
         do k=1,3
           COSA(k)=cos(LAT%ANG(k)*deg)
           SINA(k)=sin(LAT%ANG(k)*deg)
           AX(k)=LAT%AX(k)
         enddo
-        G(1,1)=(AX(2)*AX(3)*SINA(1))**2/VOL**2
-        G(2,2)=(AX(1)*AX(3)*SINA(2))**2/VOL**2
-        G(3,3)=(AX(1)*AX(2)*SINA(3))**2/VOL**2
-        G(1,2)=AX(1)*AX(2)*AX(3)**2*(COSA(1)*COSA(2)-COSA(3))/VOL**2
+        G(1,1)=(AX(2)*AX(3)*SINA(1))**2/VOL2
+        G(2,2)=(AX(1)*AX(3)*SINA(2))**2/VOL2
+        G(3,3)=(AX(1)*AX(2)*SINA(3))**2/VOL2
+        G(1,2)=AX(1)*AX(2)*AX(3)**2*(COSA(1)*COSA(2)-COSA(3))/VOL2
         G(2,1)=G(1,2)
-        G(1,3)=AX(1)*AX(3)*AX(2)**2*(COSA(1)*COSA(3)-COSA(2))/VOL**2
+        G(1,3)=AX(1)*AX(3)*AX(2)**2*(COSA(1)*COSA(3)-COSA(2))/VOL2
         G(3,1)=G(1,3)
-        G(2,3)=AX(2)*AX(3)*AX(1)**2*(COSA(2)*COSA(3)-COSA(1))/VOL**2
+        G(2,3)=AX(2)*AX(3)*AX(1)**2*(COSA(2)*COSA(3)-COSA(1))/VOL2
         G(3,2)=G(2,3)
       else
         G=0.D0
       endif
+	  LAT%G = G
       end subroutine GET_METRIC
 
 

src/simath.f

@@ -90,7 +90,7 @@ real(kind(1.D0)) function  DW_PHI1(x)
         dx=x/100
         DO i=2,100
           ksi=(i-1)*dx
-          Z=Z+P(MOD(i,2))*ksi/(exp(ksi)-1.D0)
+          Z=Z+P(MOD(i,2)+1)*ksi/(exp(ksi)-1.D0)
         ENDDO
         DW_PHI1=Z*dx/3.D0
       endif