hdiff output

r32758/orbitals.f90 2017-06-09 15:30:17.319986682 +0100 r32757/orbitals.f90 2017-06-09 15:30:17.539989744 +0100
  5:   5: 
  6:     SUBROUTINE ORBITALS_INIT(NORBS, NROTS)  6:     SUBROUTINE ORBITALS_INIT(NORBS, NROTS)
  7:   7: 
  8:         ! Initialise information needed for orbital optimisation.  8:         ! Initialise information needed for orbital optimisation.
  9:         ! This requires reading in the integral arrays from file.  9:         ! This requires reading in the integral arrays from file.
 10:  10: 
 11:         USE COMMONS, ONLY: R2INTS, DIPINTS, ORBVAREXPONENT 11:         USE COMMONS, ONLY: R2INTS, DIPINTS, ORBVAREXPONENT
 12:  12: 
 13:         INTEGER, INTENT(OUT) :: NORBS, NROTS 13:         INTEGER, INTENT(OUT) :: NORBS, NROTS
 14:  14: 
 15:         INTEGER :: ORB_FILE_UNIT, GETUNIT 15:         INTEGER :: ORB_FILE_UNIT
 16:  16: 
 17:         CHARACTER (LEN=10) :: WORD 17:         CHARACTER (LEN=10) :: WORD
 18:  18: 
 19:         LOGICAL :: END 19:         LOGICAL :: END
 20:  20: 
 21:         ! First check we've initialised sensibly. 21:         ! First check we've initialised sensibly.
 22:         IF (ORBVAREXPONENT.LT.1) THEN 22:         IF (ORBVAREXPONENT.LT.1) THEN
 23:             STOP 'Need to set penalty function exponent to a positive value.' 23:             STOP 'Need to set penalty function exponent to a positive value.'
 24:         END IF 24:         END IF
 25:  25: 
 26:         NORBS = -1 26:         NORBS = -1
 27:  27: 
 28:         CALL CHECK_FILE_EXISTS('INTEGRALS_INFO', EX=.TRUE.) 28:         CALL CHECK_FILE_EXISTS('INTEGRALS_INFO', EX=.TRUE.)
 29:  29: 
 30:         ORB_FILE_UNIT=GETUNIT() 30:         CALL FILE_OPEN('INTEGRALS_INFO',ORB_FILE_UNIT,.FALSE.)
 31:         OPEN(UNIT=ORB_FILE_UNIT,FILE='INTEGRALS_INFO',STATUS='OLD') 31:         CALL INPUT(END, ORB_FILE_UNIT)
 32:         READ(ORB_FILE_UNIT, *) NORBS 
 33:  32: 
  33:         CALL READI(NORBS)
 34:         NROTS = NORBS * (NORBS - 1) / 2 34:         NROTS = NORBS * (NORBS - 1) / 2
 35:         !PRINT '(A,3I10)','NORBS,NROTS,ORBVAREXPONENT=',NORBS,NROTS,ORBVAREXPONENT 35:         !PRINT '(A,3I10)','NORBS,NROTS,ORBVAREXPONENT=',NORBS,NROTS,ORBVAREXPONENT
 36:  36: 
 37:         IF (NORBS.LT.1) THEN 37:         IF (NORBS.LT.1) THEN
 38:             STOP 'Could not successfully read the system information from file.' 38:             STOP 'Could not successfully read the system information from file.'
 39:         ELSE IF (NROTS.LT.1) THEN 39:         ELSE IF (NROTS.LT.1) THEN
 40:             STOP 'We appear to have no possible orbital rotations.' 40:             STOP 'We appear to have no possible orbital rotations.'
 41:         END IF 41:         END IF
 42:  42: 
 43:         CLOSE(ORB_FILE_UNIT) 43:         CLOSE(ORB_FILE_UNIT)
 59:     SUBROUTINE READ_INTEGRALS(INT_FILENAME, NORBS, INTS) 59:     SUBROUTINE READ_INTEGRALS(INT_FILENAME, NORBS, INTS)
 60:  60: 
 61:         ! Given filename of intergral files reads the integrals within 61:         ! Given filename of intergral files reads the integrals within
 62:         ! it into the array INTS. Integrals are assumed to always have 62:         ! it into the array INTS. Integrals are assumed to always have
 63:         ! dimension NORBSxNORBS. 63:         ! dimension NORBSxNORBS.
 64:  64: 
 65:         CHARACTER (LEN=*), INTENT(IN) :: INT_FILENAME 65:         CHARACTER (LEN=*), INTENT(IN) :: INT_FILENAME
 66:         INTEGER, INTENT(IN) :: NORBS 66:         INTEGER, INTENT(IN) :: NORBS
 67:         DOUBLE PRECISION, INTENT(OUT) :: INTS(NORBS, NORBS) 67:         DOUBLE PRECISION, INTENT(OUT) :: INTS(NORBS, NORBS)
 68:         CHARACTER (LEN=32) :: CALCID 68:         CHARACTER (LEN=32) :: CALCID
 69:         INTEGER :: INT_FILE_UNIT, I, GETUNIT 69:         INTEGER :: INT_FILE_UNIT, I
 70:         LOGICAL :: END 70:         LOGICAL :: END
 71:  71: 
 72:         CALL CHECK_FILE_EXISTS(INT_FILENAME, EX=.TRUE.) 72:         CALL CHECK_FILE_EXISTS(INT_FILENAME, EX=.TRUE.)
 73:         INT_FILE_UNIT=GETUNIT() 73:         CALL FILE_OPEN(INT_FILENAME,INT_FILE_UNIT,.FALSE.)
 74:         OPEN(UNIT=INT_FILE_UNIT,FILE=INT_FILENAME,STATUS='OLD') 74:         CALL INPUT(END, INT_FILE_UNIT)
 75: !        CALL FILE_OPEN(INT_FILENAME,INT_FILE_UNIT,.FALSE.) 75:         CALL READU(CALCID)
 76: !        CALL INPUT(END, INT_FILE_UNIT) 
 77:         READ(INT_FILE_UNIT, '(A)') CALCID 
 78:         !CALL FILE_OPEN(INT_FILENAME,INT_FILE_UNIT,.FALSE.) 
 79:         !CALL INPUT(END, INT_FILE_UNIT) 
 80:         !CALL READU(CALCID) 
 81:         !PRINT '(A,I10,X,A)','INT_FILE_UNIT,CALCID=',INT_FILE_UNIT,CALCID 76:         !PRINT '(A,I10,X,A)','INT_FILE_UNIT,CALCID=',INT_FILE_UNIT,CALCID
 82:         !PRINT '(A,I10,A)','NORBS=',NORBS 77:         !PRINT '(A,I10,A)','NORBS=',NORBS
 83:  78: 
 84:         DO I = 1, NORBS 79:         DO I = 1, NORBS
 85:            READ(INT_FILE_UNIT,*) INTS(I,1:NORBS)  80:            READ(INT_FILE_UNIT,*) INTS(I,1:NORBS) 
 86:            !PRINT *,INTS(I,1:NORBS) 81:            !PRINT *,INTS(I,1:NORBS)
 87:         END DO 82:         END DO
 88:  83: 
 89:         CLOSE(INT_FILE_UNIT) 
 90:     END SUBROUTINE READ_INTEGRALS 84:     END SUBROUTINE READ_INTEGRALS
 91:  85: 
 92:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST) 86:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST)
 93:  87: 
 94:         ! Obtains the Power of the Orbital Variance locality measure for 88:         ! Obtains the Power of the Orbital Variance locality measure for
 95:         ! a given position. If GTEST is set also returns the gradient. 89:         ! a given position. If GTEST is set also returns the gradient.
 96:  90: 
 97:         ! Note that for intermediates dependent upon coefficient of current 91:         ! Note that for intermediates dependent upon coefficient of current
 98:         ! (rotated) orbital set in the original MO basis we use indexing 92:         ! (rotated) orbital set in the original MO basis we use indexing
 99:         ! of (ORIG_ORB,NEW_ORB) in all cases. 93:         ! of (ORIG_ORB,NEW_ORB) in all cases.
109:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS)103:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS)
110:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS)104:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS)
111: 105: 
112:         INTEGER :: I,J,K106:         INTEGER :: I,J,K
113: 107: 
114:         CALL GET_ROTATION(X, ROTATION)108:         CALL GET_ROTATION(X, ROTATION)
115: 109: 
116:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)110:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)
117:         DO I = 1, NORBS111:         DO I = 1, NORBS
118:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2)112:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2)
119:             !PRINT*,ORBVAR(I) 
120:             DO J = 1, NORBS113:             DO J = 1, NORBS
121:                 DO K = 1, NORBS114:                 DO K = 1, NORBS
122:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(J,I)*ROTATION(K,I)115:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(I,J)*ROTATION(I,K)
123:                 END DO116:                 END DO
124:             END DO117:             END DO
125:             !PRINT *,'ORBITAL',I,',ORBVAR=',ORBVAR(I),',DIPOLES=',ORBDIPOLES(1:3,I)118:             !PRINT *,'ORBITAL',I,',ORBVAR=',ORBVAR(I),',DIPOLES=',ORBDIPOLES(1:3,I)
126:         END DO119:         END DO
127: 120: 
128:         !PRINT *,'ORBVARS=',ORBVAR121:         !PRINT *,'ORBVARS=',ORBVAR
129:         !PRINT *,'ORBVARSEXP=',ORBVAR**ORBVAREXPONENT122:         !PRINT *,'ORBVARSEXP=',ORBVAR**ORBVAREXPONENT
130: 123: 
131:         LOCALITY = SUM(ORBVAR ** ORBVAREXPONENT)124:         LOCALITY = SUM(ORBVAR ** ORBVAREXPONENT)
132: 125: 
240: 233: 
241:         DO TRIIND = 1, NROTS234:         DO TRIIND = 1, NROTS
242:             CALL DECODE_TRIIND(TRIIND, I, J)235:             CALL DECODE_TRIIND(TRIIND, I, J)
243:             IF (TRIIND.EQ.DERIVTRIIND) THEN236:             IF (TRIIND.EQ.DERIVTRIIND) THEN
244:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)237:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)
245:             ELSE238:             ELSE
246:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)239:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)
247:             END IF240:             END IF
248:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT)241:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT)
249:         END DO242:         END DO
 243: 
250:         !DO TRIIND = 1, NORBS244:         !DO TRIIND = 1, NORBS
251:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS)245:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS)
252:         !END DO246:         !END DO
253: 247: 
254:     END SUBROUTINE GET_ROTATION_DERIVATIVE248:     END SUBROUTINE GET_ROTATION_DERIVATIVE
255: 249: 
256:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT)250:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT)
257: 251: 
258:         ! Obtains Givens rotation between orbitals I and J with angle theta.252:         ! Obtains Givens rotation between orbitals I and J with angle theta.
259: 253: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0