From e3c0accbcb94a380a5fecfe0653e004f8a2b0968 Mon Sep 17 00:00:00 2001
From: Sylvain Tricot <sylvain.tricot@univ-rennes1.fr>
Date: Wed, 9 Feb 2022 12:22:05 +0100
Subject: [PATCH] Added and updated coefpq.f

The file coefpq.f was updated to be compatible with Python
bindings. The common /Q_ARRAY/ was refactored as a module in
memalloc/modules.f and allocated in memalloc/allocation.f.
---
 src/msspec/spec/fortran/memalloc/allocation.f |  2 +
 src/msspec/spec/fortran/memalloc/modules.f    | 14 ++++++
 .../fortran/phd_ce_noso_nosp_nosym/coefpq.f   | 46 +++++++++++++++++++
 3 files changed, 62 insertions(+)
 create mode 100644 src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coefpq.f

diff --git a/src/msspec/spec/fortran/memalloc/allocation.f b/src/msspec/spec/fortran/memalloc/allocation.f
index 1ce3b3d..932aa3e 100644
--- a/src/msspec/spec/fortran/memalloc/allocation.f
+++ b/src/msspec/spec/fortran/memalloc/allocation.f
@@ -25,6 +25,7 @@
         USE OUTUNITS_MOD
         USE PARCAL_MOD
         USE PARCAL_A_MOD
+        USE Q_ARRAY_MOD
         USE RELADS_MOD
         USE RELAX_MOD
         USE RESEAU_MOD
@@ -136,6 +137,7 @@
         CALL ALLOC_OUTUNITS()
         CALL ALLOC_PARCAL()
         CALL ALLOC_PARCAL_A()
+        CALL ALLOC_Q_ARRAY()
         CALL ALLOC_RELADS()
         CALL ALLOC_RELAX()
         CALL ALLOC_RENORM()
diff --git a/src/msspec/spec/fortran/memalloc/modules.f b/src/msspec/spec/fortran/memalloc/modules.f
index 9e8ab0d..bc796d9 100644
--- a/src/msspec/spec/fortran/memalloc/modules.f
+++ b/src/msspec/spec/fortran/memalloc/modules.f
@@ -417,6 +417,20 @@ C=======================================================================
         END SUBROUTINE ALLOC_PARCAL_A
       END MODULE PARCAL_A_MOD
 
+C=======================================================================
+      MODULE Q_ARRAY_MOD
+        IMPLICIT NONE
+        REAL, ALLOCATABLE, DIMENSION(:) :: Q
+      CONTAINS
+        SUBROUTINE ALLOC_Q_ARRAY()
+          USE DIM_MOD
+          IF (ALLOCATED(Q)) THEN
+            DEALLOCATE(Q)
+          ENDIF
+          ALLOCATE(Q(NGR_M))
+        END SUBROUTINE ALLOC_Q_ARRAY
+      END MODULE Q_ARRAY_MOD
+
 C=======================================================================
       MODULE RELADS_MOD
         IMPLICIT NONE
diff --git a/src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coefpq.f b/src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coefpq.f
new file mode 100644
index 0000000..04413c8
--- /dev/null
+++ b/src/msspec/spec/fortran/phd_ce_noso_nosp_nosym/coefpq.f
@@ -0,0 +1,46 @@
+C
+C======================================================================
+C
+      SUBROUTINE COEFPQ(NAT,NGR)
+C
+C  This subroutine computes the P(n,m) and Q(n) coefficients 
+C   involved in the correlation expansion formulation
+C
+C  Reference  : equations (2.15) and (2.16) of 
+C               H. Zhao, D. Sebilleau and Z. Wu, 
+C               J. Phys.: Condens. Matter 20, 275241 (2008)
+C
+C                                                    H.-F. Zhao 2007
+C
+      USE DIM_MOD
+      USE Q_ARRAY
+C
+      INTEGER NAT,NGR
+C
+      REAL CMN(NGR_M,NGR_M),P(NGR_M,NGR_M)
+C
+C
+      IF(NGR.GT.NAT) THEN
+        WRITE(6,*) 'NGR is larger than NAT, which is wrong'
+        STOP
+      ENDIF
+C
+      CALL CMNGR(NAT,NGR,CMN)
+C
+      DO N=1,NGR
+        P(N,N)=1.
+        Q(N)=P(N,N)
+        DO M=N+1,NGR
+          P(N,M)=0.
+          DO I=N,M-1
+            P(N,M)=P(N,M)-P(N,I)*CMN(I,M)
+          ENDDO
+          Q(N)=Q(N)+P(N,M)
+C
+        ENDDO
+C
+      ENDDO
+C
+      RETURN
+C
+      END