[Wien] Problem with mixer

Björn Hülsen huelsen at fhi-berlin.mpg.de
Mon Oct 22 16:07:31 CEST 2007


Dear L .Marks,

thanks for the code.
But, as I wrote in my first email, I used the default input for the mixer with a history setting of 8.

MSEC1  0.0   YES  (BROYD/PRATT, extra charge (+1 for additional e), norm)
0.40            mixing FACTOR for BROYD/PRATT scheme
1.00  1.00      PW and CLM-scaling factors
9999  8         idum, HISTORY

Best regards
Bjoern

Laurence Marks wrote:
> There is a small bug in one of the mixer files if you ask for 15 or
> more limited memory steps which is fixed in the attached file (which
> also has a couple of very minor improvements). This only occurs if you
> ask for a large number, and there is NO EVIDENCE that this is good. I
> strongly recommend that you don't use more than the recommended 6-10,
> perhaps 8 as a default.
> 
> On 10/22/07, Stefaan Cottenier <Stefaan.Cottenier at fys.kuleuven.be> wrote:
>>> But twice I encountered a problem with the mixer which aborted the SCF cycle with the following error message:
>>>
>>>
>> I can offer no solution, just want to tell that when reading your
>> description I realized I have met the same problem a few times during
>> the past weeks. Similarly to you, I ignored it, blaming the notoriously
>> unstable file system of that computer. Just restarting was indeed possible.
>>
>> Stefaan
>>
>>
>>> forrtl: severe (24): end-of-file during read, unit 32, file /scr/thlc999/huelsen/cms/37_afm_tet/37_afm_tet.broyd2001
>>>
>>> The requested file and some other *broyd* were empty:
>>>
>>>  > ll *broyd*
>>> -rw-r--r--    1 huelsen  th             24 Oct 20 02:17 37_afm_tet.broyd1
>>> -rw-r--r--    1 huelsen  th              0 Oct 19 23:04 37_afm_tet.broyd2
>>> -rw-r--r--    1 huelsen  th              0 Oct 20 01:43 37_afm_tet.broyd2001
>>> -rw-r--r--    1 huelsen  th              0 Oct 20 01:55 37_afm_tet.broyd2002
>>> -rw-r--r--    1 huelsen  th              0 Oct 20 02:06 37_afm_tet.broyd2003
>>> -rw-r--r--    1 huelsen  th              0 Oct 20 02:17 37_afm_tet.broyd2004
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 19 23:50 37_afm_tet.broyd2005
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 00:01 37_afm_tet.broyd2006
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 00:13 37_afm_tet.broyd2007
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 00:24 37_afm_tet.broyd2008
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 00:35 37_afm_tet.broyd2009
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 00:47 37_afm_tet.broyd2010
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 00:58 37_afm_tet.broyd2011
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 01:10 37_afm_tet.broyd2012
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 01:21 37_afm_tet.broyd2013
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 01:32 37_afm_tet.broyd2014
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 01:43 37_afm_tet.broyd2015
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 01:55 37_afm_tet.broyd2016
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 02:06 37_afm_tet.broyd2017
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 20 02:17 37_afm_tet.broyd2018
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 19 22:42 37_afm_tet.broyd2019
>>> -rw-r--r--    1 huelsen  th           3.5M Oct 19 22:53 37_afm_tet.broyd2020
>>>
>>> The first time this error occurred we had some troubles with the filesystem. So I ignored it and restarted the calculation that afterwards converged successfully.
>>> But the second time the error occurred during another calculation on a different compute cluster that runs fine. Other calculations run smoothly on it.
>>> So this might be a problem stemming from Wien2k. In both cases I used the default settings in case.inm.
>>> Do you have any idea what might cause this problem?
>>>
>>> Best regards
>>> Bjoern
>>>
>>>
>>>
>>
>> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
>>
>> _______________________________________________
>> Wien mailing list
>> Wien at zeus.theochem.tuwien.ac.at
>> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
>>
> 
> 
> 
> ------------------------------------------------------------------------
> 
>         SUBROUTINE qmix7(clmnew,clmold,roknew,rokold,jspin,nq3,nat,lmmax,jri,qmx, &
>                  ndm,jatom1,ll,alx,aly,alz,dmat,alx_old,aly_old,alz_old,dmat_old, &
>                  mult,inversion,scl1,scl2,nbroyd,nmax,scale,LM,ndif, KVECVL, &
>                  mode,VI,MIX, qmx_input,qtot,zz,fname )
>         IMPLICIT REAL*8 (A-H,O-Z)
>         INCLUDE 'param.inc'
> !                                                                       
>         INTEGER, INTENT(IN)       :: ndm,jatom1(nat,4,3),ll(nat,4,3)
>         REAL*8, INTENT(IN)        :: alx_old(nat,4,3),aly_old(nat,4,3),alz_old(nat,4,3)
>         COMPLEX*16, INTENT(IN)    :: dmat_old(nat,4,-3:3,-3:3,3)
>         REAL*8, INTENT(INOUT)     :: alx(nat,4,3),aly(nat,4,3),alz(nat,4,3)
>         COMPLEX*16, INTENT(INOUT) :: dmat(nat,4,-3:3,-3:3,3)
>         COMPLEX*16 ROKOLD,ROKNEW
>         INTEGER LM(2,NCOM,NAT)
>         DIMENSION CLMOLD(NRAD,NCOM,NAT,2),CLMNEW(NRAD,NCOM,NAT,2),       &
>         ROKNEW(nq3,2),ROKOLD(nq3,2), dmixout(4)
>         INTEGER mult(*)
>         DIMENSION LMMAX(NAT),JRI(NAT),SCALE(*)
>         integer  KVECVL(3,nq3)
> !
>         CHARACTER*80 fname,fname2
>         CHARACTER*5 MIX
>         real*8 ZZ(*)
>         real*8,allocatable :: PS(:,:),SB(:), Y(:,:),S(:,:)
>         real*8,allocatable :: YY(:,:),SC1(:), FCURRENT(:),PCURRENT(:)
>         real*8,allocatable :: SS(:,:)
>         real*8,allocatable :: FHIST(:), PWHIST(:), CLMHIST(:)
>         logical inversion,reset,DynScale, there
>         logical DoMSec, NoLocal
>         character*4 MODUS, VERSION
>         parameter (VERSION='VER3')
>         parameter (ACUT=138.6D0, ONE=1.D0, ZERO=0.0D0)
> !                                                                       
> !       CHARGE DENSITY MIXING USING Various Methods                        
> !       Authors, L. D. Marks and R. Luke
> !       Version 3 (Beta +), May 2007
> !
> !       Update History
> !       May 20th
> !               Increased lower limit in Pratt to 0.01 -- otherwise below Wolf condition
> !               Cleaned up a couple of format statements
> !       May 21st
> !               Lined up format 1002
> !               Added back writing of :RED (sreduction, qmx)
> !               Corrected reduction -- should not matter
> !       May 24th
> !               Gone to full internal control
> !
> !       Read debug controls
>         rtrap=0.15
>         IDSCALE=0
>         iuse=0
>         inquire(file='Tester', exist=there)
>         if(there)then
>                 open(unit=1,file='Tester')
>                 read(1,*,err=983,end=983)IDTMP,rtrap
>                 IDSCALE=IDTMP
>                 goto 984
> 983             IDSCALE=0
>                 rtrap =0.15
> 984             continue
>         endif
> !       Traps for consistency (and later)
>         iuse=0
>         DoMSec=.true.
> !
> !       Output the Options
>         write(21,2001)' '
>         write(21,2001)'******************************************************'
>         write(21,2001)'* MULTISECANT MIXING OPTIONS                         *'
>         if(DoMSec) &
>         write(21,2001)'* NFL Multiscant calculated                          *'
>         if(rtrap .gt. 0.) &
>         write(21,2001)'* NFL Limit ',rtrap,'                               *'
>         if(idscale .eq. 0) &
>         write(21,2001)'* Hybrid Charge/PW/CLM Limits Used                   *'
>         if(idscale .eq. 1) &
>         write(21,2001)'* Model Independent Scaling Limits Used              *' 
>         if(idscale .eq. 2) &
>         write(21,2001)'* Old and New Scaling Limits Used                    *' 
>         nmaxx=16
>         if(nmax.ge.nmaxx) then
>          write(21,'("* NMAX gt ",i2," not supported and automatically reduced *")') nmaxx
>          nmax=nmaxx-1
>         endif
>         write(21,2002)'* Max Number of Memory Steps ',NMAX,'                    *'
>         write(21,2001)'******************************************************'
>         write(21,2001)' '
> 2002    format(a,i4,a)
> !       Determine the total Size
> !
>         MAXMQ=NQ3
>         idmat=0
>         DO idm=1,ndm
>                 DO  JATOM=1,NAT
>                         IF(jatom1(jatom,1,idm).EQ.0) EXIT
>                         DO iorb=1,4
>                                 IF(jatom1(jatom,iorb,idm).EQ.0) EXIT      
>                                 l=ll(jatom,iorb,idm)
>                                 idmat=idmat+3+(l*l+1)*(l*l+1)*2
>                         ENDDO
>                 ENDDO
>         ENDDO
>         if(inversion)then
>                 NBSIZE=jspin*nq3+NRAD*SUM(LMMAX(1:NAT))*jspin+idmat
>         else
>                 NBSIZE=jspin*nq3*2+NRAD*SUM(LMMAX(1:NAT))*jspin+idmat
>         endif
>         allocate ( PS(NBSIZE,2))
> !
> !       Scaling for Plane waves -- has an effect !!!!
>         scl_plane=sqrt(VI)*scl1
> !
> !       Pack values into PS Array
>         call PackValues(PS,ROKNEW,ROKOLD,CLMOLD,CLMNEW, &
>              alx_old,aly_old,alz_old,alx,aly,alz,dmat,dmat_old, &
>              NPLANE,JSPIN,NAT,JRI,LMMAX,JATOM,NDM,jatom1,ll,nq3,              &
>              SPLANE, SCHARGE, scl_plane, MAXMIX,NBSIZE,MaxMQ,inversion)
> !
> !       Splane is the mean-square distance in real space
> !       The factor of ndif puts it on a comparable scale to :DIS, :ENE
>         if(.not.inversion)Splane=Splane*0.5D0
>         if(Splane.gt.1d-20)Splane=sqrt(Splane)*ndif
>         write(21,210)':PLANE:  INTERSTITIAL DISTAN ',Splane
>         write(21,210)':CHARG:  CLM CHARGE   DISTAN ',SCharge
> !
> !       Apply general bounds based upon residues and type of atoms
>         call stepbound(sreduction,nat,zz,qmx_input,qmx,qtot,splane,scharge)
>         write(21,4141)sreduction,qmx
> !
> !--------------------------------------------------------------------
> !       Start Collecting History about previous steps
> !       Arrange into S, Y matrices as we go
> !--------------------------------------------------------------------
> !
>         read(31,err=6001,end=6001)dmix,niter,modus
>         dmix_last=dmix
>         goto 6002
> 6001    niter=0
> 6002    niter=niter+1
>         MEMALL=MIN(NMAX+1,NITER)
>         ISKIP=NITER-MEMALL
>         MEMORY=MEMALL-1
>         allocate (SB(MAXMIX))
>         allocate (FCURRENT(MAXMIX), PCURRENT(MAXMIX))
> !
> !       Current point, not most efficient
>         PCURRENT(1:MAXMIX)=PS(1:MAXMIX,1)
>         FCURRENT(1:MAXMIX)=PS(1:MAXMIX,2)-PCURRENT(1:MAXMIX)
>         deallocate (PS)
> !
>         allocate (FHIST(MEMALL), PWHIST(MEMALL), CLMHIST(MEMALL))
> !
>         FHIST(MEMALL)  =dot_product(FCURRENT,FCURRENT)
>         PWHIST(MEMALL) =dot_product(FCURRENT(1:NPLANE),FCURRENT(1:NPLANE))
>         CLMHIST(MEMALL)=FHIST(MEMALL)-PWHIST(MEMALL)
>         if(MEMORY .gt. 0)then
>                 allocate (Y(MAXMIX,MEMORY), S(MAXMIX,MEMORY))
>                 allocate (YY(MEMORY,MEMORY),SS(MEMORY,MEMORY))
>                 allocate (SC1(MEMORY))
> !PB
>                 close (32)
>                 if(iskip.gt.0)write(21,*)'Skipping ',iskip,' Total ',niter
> !                DO N=1,ISKIP
> !!       Note: use SB as workspace
> !                   read(32)SB
> !                   read(32)SB
> !                enddo
>                 DO N=1,MEMORY
>                    call open_broydfile(fname,iskip+n,fname2)
> !                   write(21,*) 'open ',iskip,n,fname2
>                    open(32,file=fname2,status='old',form='unformatted')
>                    read(32)SB
>                    S(1:MAXMIX,N)=PCURRENT(1:MAXMIX)-SB(1:MAXMIX)
>                    read(32)SB
>                    close (32)
>                    Y(1:MAXMIX,N)=FCURRENT(1:MAXMIX)-SB(1:MAXMIX)
>                    FHIST(N)  =dot_product(SB(1:MAXMIX),SB(1:MAXMIX))
>                    PWHIST(N) =dot_product(SB(1:NPLANE),SB(1:NPLANE))
>                    CLMHIST(N)=FHIST(N)-PWHIST(N)
>                 enddo
>         endif
> !
>                    call open_broydfile(fname,iskip+memory+1,fname2)
> !                   write(21,*) 'open ',iskip,memory,fname2
>                    open(32,file=fname2,status='unknown',form='unformatted')
>         write(32)PCURRENT
>         write(32)FCURRENT
>         close(32)
>                    if(niter-nmaxx.ge.0) then
>                       call open_broydfile(fname,niter-nmaxx+1,fname2)
>                       OPEN (32,FILE=FNAME2,STATUS='REPLACE')
>                       close (32,STATUS='DELETE')
> !     print*, 'remove ',fname2, niter
>                    endif
> !
> !--------------------------------------------------------------------
> !       End Collecting History about previous steps
> !--------------------------------------------------------------------
> !
>         if(MEMORY .eq. 0)goto 119
> !
> !--------------------------------------------------------------------
> !       Renormalize the S & Y vectors
> !--------------------------------------------------------------------
> !
>         do J=1,MEMORY
>            T1=0.
> !
>            DO K=1,MAXMIX
>                 T1=T1+Y(K,J)*Y(K,J)
>            enddo
> !       Renormalize
>            T1=1.D0/sqrt(T1)
>            DO K=1,MAXMIX
>                 S(K,J)=S(K,J)*T1
>                 Y(K,J)=Y(K,J)*T1
>            ENDDO
> !
>         enddo
> !
> !--------------------------------------------------------------------
> !       End Renormalizing the S & Y vectors
> !--------------------------------------------------------------------
> !
> !--------------------------------------------------------------------
> !       Generate the MEMORY x MEMORY Matrices
> !       Also generate scaling information
> !--------------------------------------------------------------------
> !
> !       Note: should use dgemms here -- for later
> !
>         YYSUM1=0.0
>         DO J=1,MEMORY
>           DO K=1,J
>                 SS(J,K) =dot_product(S (1:MAXMIX,J),S (1:MAXMIX,K))
>                 YY(J,K) =dot_product(Y (1:MAXMIX,J),Y (1:MAXMIX,K))
>           enddo
> !
> !       Some scaling type material, not yet completely used
> !       These measure how well previous steps have performed
> !       1D-30 term added to avoid divisions by zero
> !
>           SC1(J)=sqrt(SS(J,J)*YY(J,J)/(YY(J,J)*YY(J,J)+1D-30))
> 
>           YYSUM1=YYSUM1+YY (J,J)
> !       Running averages
>           if(J .eq. 1)then
>                 ASCL1=SC1(J)
>           else
>                 ASCL1=(ASCL1*2+SC1(J))/3.D0
>           endif
>         enddo
> !       Do transpose part as well
>         DO J=1,MEMORY
>                 DO K=1,J-1
>                         SS(K,J)=SS(J,K)
>                         YY(K,J)=YY(J,K)
>                 enddo
>         enddo
> !
> !       Get the Condition Number of the matrices
> !       Useful diagnostic, although these are horribly big!
>         DIAG=1D-3
>         CondYY =GetCondNo(YY ,MEMORY,DIAG,iYY)
>         CondSS =GetCondNo(SS ,MEMORY,DIAG,ISS)
> !
> 
>         write(21,1002) ' YY,SS Cond No',CondYY,CondSS       
> !
> !       Look at the seperate Plane and RMT progress scalings
>         YPL=dot_product(Y(1:NPLANE,MEMORY),Y(1:NPLANE,MEMORY))
>         SPL=dot_product(S(1:NPLANE,MEMORY),S(1:NPLANE,MEMORY))
>         SCP=sqrt(SPL*YPL/(YPL*YPL+1D-50))
>         YRM=YY(MEMORY,MEMORY)-YPL
>         SRM=SS(MEMORY,MEMORY)-SPL
>         SCR=sqrt(SRM*YRM/(YRM*YRM+1D-50))
> !
> !
> !--------------------------------------------------------------------
> !       Bound the step here
>         call LimitDMIX(Y,YY,FCURRENT,FHIST,PWHIST,CLMHIST,MAXMIX,MEMALL,MEMORY, &
>                              ascl1, qmx_input, dmix_last, qmx, dmixout, &
>                              IDSCALE,Nplane, PM1, rtrap )
> !  
>         DMIXM=dmixout(1)
>         DMIX=DMIXM
>         qmx=dmixm
> !
> !--------------------------------------------------------------------
> !  
>         write(21,1002)' SS/YY Scales ',SC1(MEMORY),ASCL1,SCP,SCR
> !
> !--------------------------------------------------------------------
> !                                                                       
> !       Is it the current version or not?
>         if(MODUS .ne. VERSION) goto 119
> !
> !--------------------------------------------------------------------
> !       Calculate the New Steps, currently with 3 modes             
> !
> !--------------------------------------------------------------------
> !       Set Diagonal Initial Matrix
> !        call SetHZero(HZERO,DMIX,SCL1,MAXMIX,NPLANE)
> !--------------------------------------------------------------------
> !
> !      Note: Y and S referenced to current position here -- works better
>        call  MSEC1(Y,S,YY,FCURRENT,SB,MAXMIX,MEMORY,DMIXM,IFAIL)
> !      Inversion bombed, go to Pratt
>        if(IFAIL .ne. 0)goto 119
>        write(21,1002)' DMIXM and Projections   ',DMIXM,PM1
> !
> !       Which one to use -- default MSEC1
>         MIX='MSEC1'
>         SB(1:MAXMIX)=PCURRENT(1:MAXMIX)+SB(1:MAXMIX)
>         DMIXUSED=DMIXM
> !
>         DMIX=DMIXUSED
>         QMX=DMIX
> !
> !--------------------------------------------------------------------
> !       Done with calculating the updates
> !
> !       Consistency check
>         nreduce=0
>         reduction=1.D0
> 120     DOT=0
>         BMOD=0
>         PMOD=0
> !
> !       Calculate the angle between the Pratt & Broyden Steps
> !       Also calculate the relative sizes
>         DO K=1,NPLANE
>                 SPRATT=FCURRENT(K)*DMIX 
>                 SBROY=(SB(K)-PCURRENT(K)) 
>                 DOT=SPRATT*SBROY+DOT
>                 BMOD=BMOD+SBROY*SBROY
>                 PMOD=PMOD+SPRATT*SPRATT
>         enddo
> !       Adjust to put plane waves on appropriate scale
>         DOTP=DOT 
>         BMODP=BMOD 
>         PMODP=PMOD 
>         DTTP=DOT/sqrt(BMOD*PMOD)
>         if(DTTP .gt. 1.0D0)then
>                 DTTP=0.0D0
>         else if(DTTP .lt. -1.0D0) then
>                 DTTP=180.D0
>         else
>                 DTTP=acos(DTTP)*90./asin(1.D0)
>         endif
>         DO K=1+NPLANE,MAXMIX
>                 SPRATT=FCURRENT(K)*DMIX 
>                 SBROY=(SB(K)-PCURRENT(K)) 
>                 DOT=SPRATT*SBROY+DOT
>                 BMOD=BMOD+SBROY*SBROY
>                 PMOD=PMOD+SPRATT*SPRATT
>         ENDDO
> !       Traps based upon the relative sizes and the angle
> !
> !       Test the angle
> !       If it is too big Broyden method may have collapsed
>         BMOD=sqrt(BMOD)
>         PMOD=sqrt(PMOD)
>         DTT=DOT/(BMOD*PMOD)
>         if(DTT .gt. 1.0D0)then
>                 DTT=0.0D0
>         else if(DTT .lt. -1.0D0) then
>                 DTT=180.D0
>         else
>                 DTT=acos(DTT)*90./asin(1.D0)
>         endif
> !
>         IF(abs(DTT).GT.ACUT.and.niter.ge.4)Then
>                         Write(21,98)BMOD,PMOD,DTT
>                         QMX=DMIX*MIN(2.D0*BMOD/PMOD,1.D0)
>                         QMX=min(QMX,0.075D0)*exp(-1.5*SPLANE)
>                         Goto 119
>         endif
> 
> !       Too large a step; use information about fraction projected onto Prior (PROJL)
> !       to help control this.
>         SBOUND=EXP(-2.d0*PM1)
>         BOUND1=SBOUND*PMOD/dmix                  ! 0.5 of the full step
>         BOUND2=PMOD*10.D0*SBOUND                 ! 5 times the current reduced Pratt step
> !
> !       First Bound, the smaller of these
>         BLIMIT=min(BOUND1,BOUND2)
> !
> !       If the user has a large step, be reasonable about trying it
>         BLIMIT=max(BLIMIT, PMOD*3.D0)
> !
> !       Ensure that we don't go too small
>         BLIMIT=max(BLIMIT,2.5D-2)
> !
> !       Skip trap if broyden is very small
>         if(bmod .lt. 1.D-4)goto 122
> !
>         if(BMOD .gt. BLIMIT)then
> !                       Reduce the step
>                         nreduce=nreduce+1
> !                       In case something went wrong
>                         if(nreduce .gt. 1)then
>                                 X1=0.95D0
>                         else
>                                 X1=0.99999999D0*BLIMIT/BMOD
>                         endif
>                         X2=(1.D0-X1)
>                         reduction=reduction*X1
>                         do K=1,maxmix
>                                 SB(K)=SB(K)*X1+PCURRENT(K)*X2
>                         enddo
>                         goto 120
>         endif
> !
> !       Done with traps, move to convert values back
> 122     continue
>         write(21,982)sqrt(Bmodp),sqrt(pmodp),dttp
>         if(reduction .lt. 0.9)then
>             Write(21,981)BMOD,PMOD,DTT,reduction
>         else
>             Write(21,98) BMOD,PMOD,DTT
>         endif
> 
> !--------------------------------------------------------------------
> !       Put positions back into the arrays that they came from
> !--------------------------------------------------------------------
> !
> 100     CONTINUE
>         Call UnPackValues(ROKNEW,CLMNEW, alx,aly,alz,dmat, &
>              NPLANE,JSPIN,NAT,JRI,LMMAX,JATOM,NDM,jatom1,ll,nq3,              &
>              MAXMIX,SB, scl_plane,MaxMQ,inversion)
>         rewind 31
>         write(31)DMIX,NITER,VERSION                                                         
> 
>         GOTO 121
> !                                                          
> !--------------------------------------------------------------------
> !       PRATT MIXING FOR FIRST ITERATION, after angle gets too big
> !       or something else goes wrong
> !--------------------------------------------------------------------
> !                                  
> 119     CONTINUE
> !       Note: We are resetting here, but this might not always be right
> !       If we get here via a trap on the angle perhaps we should not reset
> !
>         niter=1
>         DMIX=min(QMX,0.075D0*exp(-1.5*SPLANE))
>         DMIX=max(DMIX,0.01D0)
>         MIX='PRATT'
>         QMX=DMIX
>         DMIXUSED=DMIX
> !
>         DO K=1,MAXMIX
>                 SB(K)=PCURRENT(K)+DMIX*FCURRENT(K)
>         ENDDO
> !       Jump back to main output code
>         GOTO 100                                     
> 
> 121     CONTINUE                                                          
> !       Set QMX for output
>         DMIX=DMIXUSED
>         QMX =DMIX
> !
> !       Cleanup
>         deallocate (FCURRENT, PCURRENT,SB)
>         deallocate (FHIST, PWHIST, CLMHIST)
>         if(MEMORY .gt. 0)deallocate (Y,S,YY,SS,SC1)
> !
>         RETURN                                                            
> !                                                                       
> 200     format(':PRATT STEP USED')
> 210 	FORMAT(A,F11.7)
> 982     format(':DIRP :  |BROYD|=',D10.3,' |PRATT|=',D10.3, &
>         ' ANGLE=',f6.1,' DEGREES')
> 981     format(':DIRB :  |BROYD|=',D10.3,' |PRATT|=',D10.3, &
>         ' ANGLE=',f6.1,' DEGREES, STEP=',D10.3)
> 98      format(':DIRB :  |BROYD|=',D10.3,' |PRATT|=',D10.3, &
>         ' ANGLE=',f6.1,' DEGREES')
> 1002    format(':DIRE : ',a,10D11.3)
> 1102    format(':DIRE : ',a,4I4,' out of ',i4)
> 2001    format(a,F10.6,a)
> 2100    format(a,F10.6,2D14.6)
> 1022    FORMAT(10X,'DENSITY FOR ITERATION',I4,' PREPARED')
> 4141    format(':REDuction and DMIX in Broyd:',3f10.4,E14.5)
>         END
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Wien mailing list
> Wien at zeus.theochem.tuwien.ac.at
> http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien

-- 
Björn Hülsen
Fritz-Haber-Institute
Theory Department
Faradayweg 4-6
14195 Berlin

phone:  +49 30 8413-4863
email:  huelsen at fhi-berlin.mpg.de
www:    www.fhi-berlin.mpg.de


More information about the Wien mailing list