Line data    Source code

       1             : MODULE m_a_pulay
       2             :   USE m_juDFT
       3             :   PRIVATE
       4             :   REAL :: distance(3)
       5             :   PUBLIC a_pulay
       6             : CONTAINS
       7           0 :   SUBROUTINE a_pulay(alpha,fm,sm,l_dfpt)
       8             :     USE m_pulay
       9             :     USE m_types_mat
      10             :     USE m_types_mixvector
      11             :     USE m_mixing_history
      12             :     IMPLICIT NONE
      13             :     REAL,INTENT(IN)                 :: alpha
      14             :     TYPE(t_mixvector),INTENT(IN)    :: fm(:)
      15             :     TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
      16             :     LOGICAL,            INTENT(IN) :: l_dfpt
      17             : 
      18             :     REAL    :: fac(2)=(/1.0,0.5/)
      19             :     INTEGER,parameter :: parallel_images=2
      20             :     INTEGER,PARAMETER :: maxiter=4,mode=2
      21             :     INTEGER :: hlen,image,i,ii,local_length
      22             :     INTEGER :: local_hist(maxiter*parallel_images)
      23             : 
      24           0 :     hlen=SIZE(sm) !Number of results in history
      25           0 :     image=MOD(hlen-1,parallel_images)+1
      26           0 :     PRINT *,"HI:",hlen,image
      27           0 :     IF (hlen<parallel_images) THEN
      28             :        !create inital parallel images by simple mixing of first result
      29           0 :        sm(hlen)=sm(1)+(alpha*fac(image))*fm(1)
      30           0 :        RETURN
      31             :     ENDIF
      32             : 
      33             :     SELECT CASE(mode)
      34             :     CASE(1)
      35             :        local_length=0
      36             :        DO i=image,hlen,parallel_images
      37             :           local_length=local_length+1
      38             :           local_hist(local_length)=i
      39             :        ENDDO
      40             : 
      41             :        IF (local_length>maxiter) THEN !recreate to enable cross-image mixing
      42             :           local_length=hlen+1-parallel_images
      43             :           local_hist(1)=image
      44             :           local_hist(2:local_length)=(/(i,i=parallel_images+1,hlen)/)
      45             :        END IF
      46             :        sm(hlen)=simple_pulay(alpha*fac(image),fm(local_hist(:local_length)),sm(local_hist(:local_length)),l_dfpt)
      47           0 :        IF (hlen==parallel_images*(maxiter+1)) CALL mixing_history_limit(parallel_images)
      48             :     case(2)
      49           0 :        local_length=hlen-image
      50           0 :        local_hist(1:local_length)=(/(i,i=1,local_length)/)
      51           0 :        local_length=local_length+1
      52           0 :        local_hist(local_length)=hlen
      53           0 :        sm(hlen)=simple_pulay(alpha*fac(image),fm(local_hist(:local_length)),sm(local_hist(:local_length)),l_dfpt)
      54           0 :        IF (hlen==parallel_images*(maxiter)) CALL mixing_history_limit(parallel_images)
      55             :     END SELECT
      56             : 
      57           0 :     PRINT *,"P:",local_hist(:local_length)
      58             : 
      59             : 
      60             :   END SUBROUTINE a_pulay
      61             : 
      62           0 :   FUNCTION simple_pulay(alpha,fm,sm,l_dfpt)RESULT(sm_out)
      63             :     USE m_types_mixvector
      64             :     IMPLICIT NONE
      65             :     REAL,INTENT(IN)                 :: alpha
      66             :     TYPE(t_mixvector),INTENT(IN)    :: fm(:)
      67             :     TYPE(t_mixvector),INTENT(IN)    :: sm(:)
      68             :     LOGICAL,            INTENT(IN) :: l_dfpt
      69             :     TYPE(t_mixvector)               :: sm_out
      70             : 
      71             : 
      72           0 :     TYPE(t_mixvector)   :: mdf
      73           0 :     REAL,ALLOCATABLE    :: a(:,:),b(:),work(:)
      74           0 :     INTEGER,ALLOCATABLE :: ipiv(:)
      75             :     INTEGER             :: n,nn,info,lwork,hlen
      76             : 
      77             : 
      78             : 
      79           0 :     hlen=SIZE(sm)
      80           0 :     IF (hlen==1) THEN
      81           0 :        sm_out=sm(1)+alpha*fm(1)
      82           0 :        RETURN
      83             :     ENDIF
      84             : 
      85           0 :     ALLOCATE(a(hlen+1,hlen+1),b(hlen+1),ipiv(hlen+1),work((hlen+1)**2))
      86           0 :     DO n=1,hlen
      87           0 :        mdf=fm(n)%apply_metric(l_dfpt)
      88           0 :        DO nn=1,n
      89           0 :           a(n,nn)=mdf.dot.fm(nn)
      90           0 :           a(nn,n)=a(n,nn)
      91             :        ENDDO
      92             :     ENDDO
      93             : 
      94           0 :     a(:,hlen+1)=1.0
      95           0 :     a(hlen+1,:)=1.0
      96           0 :     a(hlen+1,hlen+1)=0.0
      97           0 :     b=0.0;b(hlen+1)=1.
      98           0 :     CALL  DSYSV( 'Upper', hlen+1, 1, a, SIZE(a,1), ipiv, b, SIZE(b,1), work, SIZE(work), INFO )
      99           0 :     sm_out=0.0*mdf
     100           0 :     DO n=1,hlen
     101           0 :        sm_out=sm_out+b(n)*(sm(n)+alpha*fm(n))
     102             :     ENDDO
     103             :   
     104           0 :   END FUNCTION SIMPLE_PULAY
     105             : 
     106             : END MODULE m_a_pulay