First attempt to Pulay...

main/mix.F90

@@ -28,6 +28,7 @@ contains
     use m_types
     use m_umix
     USE m_kerker
+    use m_pulay
     use m_types_mixvector
     USE m_distance
     use m_mixing_history
@@ -117,6 +118,8 @@ contains
     if(input%imix==0.or.it==1) CALL stmix(atoms,input,noco,fsm(it),fsm_mag,sm(it))
     !if(it>1.and.input%imix==9) CALL pulay(input%alpha,fsm,sm)
     if(it>1.and.(input%imix==3.or.input%imix==5.or.input%imix==7)) Call broyden(input%alpha,fsm,sm)
+    !PRINT *,"ATTENTION Broyden replaced by Pulay"
+    IF (it>1.and.input%imix==9) CALL pulay(input%alpha,fsm,sm)
 
     !extracte mixed density 
     inDen%pw=0.0;inDen%mt=0.0

mix/CMakeLists.txt

@@ -7,6 +7,7 @@ mix/kerker.F90
 mix/broyden_history.F90
 mix/mixing_history.F90
 mix/distance.F90
+mix/pulay.F90
 #mix/broyden2.F90
 #mix/brysh1.f90
 #mix/brysh2.f90

mix/mixing_history.F90

@@ -82,7 +82,7 @@ CONTAINS
     endif
     IF (iter_stored+1==maxiter.AND.imix.NE.9) iter_stored=0 !This is a broyden method which has to 
                                                             !be reset as soon as maxiter is reached
-    it=MIN(iter_stored+1,maxiter)
+    it=MIN(iter_stored+1,maxiter+1)
     allocate(sm(it),fsm(it))
     CALL sm(it)%alloc()
     CALL fsm(it)%alloc()
@@ -109,6 +109,7 @@ CONTAINS
     IMPLICIT NONE
     TYPE(t_mpi),INTENT(in)::mpi
     iter_stored=0
+    PRINT *, "Reset of history"
     IF (mpi%irank==0) CALL system('rm mixing_history*')
   END SUBROUTINE mixing_history_reset
   

mix/pulay.F90

+MODULE m_pulay
+  USE m_juDFT
+
+CONTAINS
+  SUBROUTINE pulay(alpha,fm,sm)
+    USE m_types_mat
+    USE m_types_mixvector
+    IMPLICIT NONE
+    REAL,INTENT(IN)                 :: alpha
+    TYPE(t_mixvector),INTENT(IN)    :: fm(:)
+    TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
+
+    ! Locals
+    INTEGER           :: h_len,n,nn
+    REAL,ALLOCATABLE  :: b(:)
+    TYPE(t_mat)       :: a
+    TYPE(t_mixvector),ALLOCATABLE :: df(:),ds(:)
+
+    h_len=SIZE(fm)-1
+    
+    sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
+    
+    PRINT *,"len:",h_len
+    DO n=1,h_len
+       PRINT *,"D:",sm(n).dot.sm(n),fm(n).dot.fm(n)
+    ENDDO
+    IF (h_len==1) RETURN
+    
+    CALL a%alloc(.TRUE.,h_len,h_len)
+
+    ALLOCATE(df(h_len),ds(h_len),b(h_len))
+
+    DO n=1,h_len
+       df(n)=fm(n+1)-fm(n)
+       ds(n)=sm(n+1)-sm(n)
+       b(n)=df(n).dot.fm(h_len+1)
+       DO nn=1,n
+          a%data_r(n,nn)=df(n).dot.df(nn)
+          a%data_r(nn,n)=a%data_r(n,nn)
+       ENDDO
+    ENDDO
+    
+    call a%inverse()
+
+    PRINT *,"b0:",b
+    b=MATMUL(a%data_r,b)
+    
+    PRINT *,"b1:",b
+    
+    DO n=1,h_len
+       sm(h_len+1)=sm(h_len+1)-b(n)*(ds(n)+alpha*df(n))
+    ENDDO
+    
+
+ 
+  END SUBROUTINE pulay
+END MODULE m_pulay