The LJ driver now does not assume upper-triangular cell

Useful for testing the rotational averaging
lab-cosmo · Sep 26, 2024 · 940de13 · 940de13
1 parent 7c17f13
commit 940de13
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 11 deletions.
diff --git a/drivers/f90/LJ.f90 b/drivers/f90/LJ.f90
@@ -189,8 +189,7 @@ SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list,
                      forces(n_list(j),:) = forces(n_list(j),:) - fij
                      pot = pot + pot_ij
                      DO k = 1, 3
-                        DO l = k, 3
-                           ! Only the upper triangular elements calculated.
+                        DO l = 1, 3
                            virial(k,l) = virial(k,l) + fij(k)*rij(l)
                         ENDDO
                      ENDDO
@@ -199,8 +198,8 @@ SUBROUTINE LJ_getall(rc, sigma, eps, natoms, atoms, cell_h, cell_ih, index_list,
                start = index_list(i) + 1
             ENDDO
 
-            ! Assuming an upper-triangular vector matrix for the simulation box.
-            volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
+            ! Works with a generic cell, even if usually it'll be upper-triuangular
+            volume = CELL_VOLUME(cell_h) !cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
             CALL LJ_longrange(rc, sigma, eps, natoms, volume, pot_lr, vir_lr)
             pot = pot + pot_lr
             DO k = 1, 3
@@ -272,8 +271,7 @@ SUBROUTINE LJMix_getall(n_type2, rc, sigma, eps, natoms, atoms, cell_h, cell_ih,
                      forces(n_list(j),:) = forces(n_list(j),:) - fij
                      pot = pot + pot_ij
                      DO k = 1, 3
-                        DO l = k, 3
-                           ! Only the upper triangular elements calculated.
+                        DO l = 1, 3
                            virial(k,l) = virial(k,l) + fij(k)*rij(l)
                         ENDDO
                      ENDDO
@@ -282,8 +280,8 @@ SUBROUTINE LJMix_getall(n_type2, rc, sigma, eps, natoms, atoms, cell_h, cell_ih,
                start = index_list(i) + 1
             ENDDO
 
-            ! Assuming an upper-triangular vector matrix for the simulation box.
-            volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
+            ! Works with a generic cell, even if usually it'll be upper-triuangular
+            volume = CELL_VOLUME(cell_h) ! cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
             CALL LJ_longrange(rc, sigma*(1-n_type2*0.4/natoms), eps*(1-n_type2*0.4/natoms), natoms, volume, pot_lr, vir_lr)
             pot = pot + pot_lr
             DO k = 1, 3

diff --git a/drivers/f90/distance.f90 b/drivers/f90/distance.f90
@@ -31,6 +31,17 @@ MODULE DISTANCE
 
       CONTAINS
 
+         DOUBLE PRECISION function CELL_VOLUME(m)
+
+         DOUBLE PRECISION m(3, 3)
+
+         CELL_VOLUME = DABS(m(1, 1) * (m(2, 2) * m(3, 3) - m(3, 2) * m(2, 3)) - &
+            m(1, 2) * (m(2, 1) * m(3, 3) - m(3, 1) * m(2, 3)) +  &
+            m(1, 3) * (m(2, 1) * m(3, 2) - m(3, 1) * m(2, 2)) )
+
+         return
+         end
+
          SUBROUTINE vector_separation(cell_h, cell_ih, ri, rj, rij, r2)
             ! Calculates the vector separating two atoms.
             !

diff --git a/drivers/f90/driver.f90 b/drivers/f90/driver.f90
@@ -32,6 +32,7 @@ PROGRAM DRIVER
          USE SG
          USE PSWATER
          USE F90SOCKETS, ONLY : open_socket, writebuffer, readbuffer, f_sleep
+         USE DISTANCE, only: CELL_VOLUME
       IMPLICIT NONE
 
       ! SOCKET VARIABLES
@@ -387,7 +388,7 @@ PROGRAM DRIVER
             enddo
          ENDIF
          isinit = .true.
-      ELSEIF (vstyle == 1) THEN
+      ELSEIF (1 == vstyle) THEN
          IF (par_count /= 3) THEN
             WRITE(*,*) "Error: parameters not initialized correctly."
             WRITE(*,*) "For LJ potential use -o sigma,epsilon,cutoff "
@@ -521,8 +522,8 @@ PROGRAM DRIVER
             ! units and storage mode used in the driver.
             cell_h = transpose(cell_h)
             cell_ih = transpose(cell_ih)
-            ! We assume an upper triangular cell-vector matrix
-            volume = cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
+            ! We compute for a generic cell, just in case (even though usually i-PI passes an upper triangular cell-vector matrix)
+            volume = CELL_VOLUME(cell_h) !cell_h(1,1)*cell_h(2,2)*cell_h(3,3)
 
             CALL readbuffer(socket, cbuf)       ! The number of atoms in the cell
             IF (verbose > 1) WRITE(*,*) "    !read!=> cbuf: ", cbuf