pulay.F90 1.51 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4
MODULE m_pulay
  USE m_juDFT

CONTAINS
5
  SUBROUTINE pulay(alpha,fm,sm,simple_steps)
Daniel Wortmann's avatar
Daniel Wortmann committed
6 7
    USE m_types_mat
    USE m_types_mixvector
8
    USE m_mixing_history
Daniel Wortmann's avatar
Daniel Wortmann committed
9 10 11
    IMPLICIT NONE
    TYPE(t_mixvector),INTENT(IN)    :: fm(:)
    TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
12 13
    INTEGER,INTENT(IN)              :: simple_steps
    REAL,INTENT(IN)                 :: alpha
Daniel Wortmann's avatar
Daniel Wortmann committed
14 15 16 17 18
    ! Locals
    INTEGER           :: h_len,n,nn
    REAL,ALLOCATABLE  :: b(:)
    TYPE(t_mat)       :: a
    TYPE(t_mixvector),ALLOCATABLE :: df(:),ds(:)
19 20
    TYPE(t_mixvector) ::mdf

Daniel Wortmann's avatar
Daniel Wortmann committed
21 22 23

    h_len=SIZE(fm)-1
    
24 25 26 27 28 29 30 31 32 33 34
    IF (h_len==0) THEN
       sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
       RETURN !No history present
    ENDIF
    IF (simple_steps>0) THEN
       IF (MOD(h_len,simple_steps).NE.0) THEN
          !Simple mixing step of periodic Pulay
          sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) 
          RETURN 
       ENDIF
    ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
35 36 37 38 39 40 41
    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)
42 43
       mdf=df(n).apply_metric()
       b(n)=mdf.dot.fm(h_len+1)
Daniel Wortmann's avatar
Daniel Wortmann committed
44
       DO nn=1,n
45
          a%data_r(n,nn)=mdf.dot.df(nn)
Daniel Wortmann's avatar
Daniel Wortmann committed
46 47 48 49 50 51 52
          a%data_r(nn,n)=a%data_r(n,nn)
       ENDDO
    ENDDO
    
    call a%inverse()

    b=MATMUL(a%data_r,b)
53
    sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
Daniel Wortmann's avatar
Daniel Wortmann committed
54 55 56
    DO n=1,h_len
       sm(h_len+1)=sm(h_len+1)-b(n)*(ds(n)+alpha*df(n))
    ENDDO
57 58
   
    IF (simple_steps>0) CALL mixing_history_limit(0)
Daniel Wortmann's avatar
Daniel Wortmann committed
59 60
  END SUBROUTINE pulay
END MODULE m_pulay