Line data    Source code

       1             : MODULE m_pulay
       2             :   USE m_juDFT
       3             : 
       4             : CONTAINS
       5           0 :   SUBROUTINE pulay(alpha,fm,sm,simple_steps,l_dfpt)
       6             :     USE m_types_mat
       7             :     USE m_types_mixvector
       8             :     USE m_mixing_history
       9             :     IMPLICIT NONE
      10             :     TYPE(t_mixvector),INTENT(IN)    :: fm(:)
      11             :     TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
      12             :     INTEGER,INTENT(IN)              :: simple_steps
      13             :     REAL,INTENT(IN)                 :: alpha
      14             :     LOGICAL,            INTENT(IN) :: l_dfpt
      15             :     ! Locals
      16             :     INTEGER           :: h_len,n,nn
      17           0 :     REAL,ALLOCATABLE  :: b(:)
      18           0 :     TYPE(t_mat)       :: a
      19           0 :     TYPE(t_mixvector),ALLOCATABLE :: df(:),ds(:)
      20           0 :     TYPE(t_mixvector) ::mdf
      21             : 
      22             : 
      23           0 :     h_len=SIZE(fm)-1
      24             : 
      25           0 :     IF (h_len==0) THEN
      26           0 :        sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
      27           0 :        RETURN !No history present
      28             :     ENDIF
      29           0 :     IF (simple_steps>0) THEN
      30           0 :        IF (MOD(h_len,simple_steps).NE.0) THEN
      31             :           !Simple mixing step of periodic Pulay
      32           0 :           sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1)
      33           0 :           RETURN
      34             :        ENDIF
      35             :     ENDIF
      36           0 :     CALL a%alloc(.TRUE.,h_len,h_len)
      37             : 
      38           0 :     ALLOCATE(df(h_len),ds(h_len),b(h_len))
      39             : 
      40           0 :     DO n=1,h_len
      41           0 :        df(n)=fm(n+1)-fm(n)
      42           0 :        ds(n)=sm(n+1)-sm(n)
      43           0 :        mdf=df(n)%apply_metric(l_dfpt)
      44           0 :        b(n)=mdf.dot.fm(h_len+1)
      45           0 :        DO nn=1,n
      46           0 :           a%data_r(n,nn)=mdf.dot.df(nn)
      47           0 :           a%data_r(nn,n)=a%data_r(n,nn)
      48             :        ENDDO
      49             :     ENDDO
      50             : 
      51           0 :     call a%inverse()
      52             : 
      53           0 :     b=MATMUL(a%data_r,b)
      54           0 :     sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
      55           0 :     DO n=1,h_len
      56           0 :        sm(h_len+1)=sm(h_len+1)-b(n)*(ds(n)+alpha*df(n))
      57             :     ENDDO
      58             : 
      59           0 :     IF (simple_steps>0) CALL mixing_history_limit(0)
      60           0 :   END SUBROUTINE pulay
      61           0 : END MODULE m_pulay