hdiff output

r32748/orbitals.f90 2017-06-08 12:30:15.112174475 +0100 r32747/orbitals.f90 2017-06-08 12:30:15.388178146 +0100
 81:            !PRINT *,INTS(I,1:NORBS) 81:            !PRINT *,INTS(I,1:NORBS)
 82:         END DO 82:         END DO
 83:  83: 
 84:     END SUBROUTINE READ_INTEGRALS 84:     END SUBROUTINE READ_INTEGRALS
 85:  85: 
 86:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST) 86:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST)
 87:  87: 
 88:         ! Obtains the Power of the Orbital Variance locality measure for 88:         ! Obtains the Power of the Orbital Variance locality measure for
 89:         ! a given position. If GTEST is set also returns the gradient. 89:         ! a given position. If GTEST is set also returns the gradient.
 90:  90: 
 91:         ! Note that for intermediates dependent upon coefficient of current 
 92:         ! (rotated) orbital set in the original MO basis we use indexing 
 93:         ! of (ORIG_ORB,NEW_ORB) in all cases. 
 94:  
 95:         USE COMMONS, ONLY: NROTS, NORBS, R2INTS, DIPINTS, ORBVAREXPONENT 91:         USE COMMONS, ONLY: NROTS, NORBS, R2INTS, DIPINTS, ORBVAREXPONENT
 96:  92: 
 97:         DOUBLE PRECISION, INTENT(IN) :: X(NROTS) 93:         DOUBLE PRECISION, INTENT(IN) :: X(NROTS)
 98:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(NROTS), LOCALITY 94:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(NROTS), LOCALITY
 99:         LOGICAL, INTENT(IN) :: GTEST 95:         LOGICAL, INTENT(IN) :: GTEST
100:  96: 
101:         DOUBLE PRECISION :: ROTATION(NORBS,NORBS), ORBVAR(NORBS) 97:         DOUBLE PRECISION :: ROTATION(NORBS,NORBS), ORBVAR(NORBS)
102:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 98:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS)
103:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS) 99:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS)
104:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS)100:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS)
109: 105: 
110:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)106:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)
111:         DO I = 1, NORBS107:         DO I = 1, NORBS
112:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2)108:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2)
113:             DO J = 1, NORBS109:             DO J = 1, NORBS
114:                 DO K = 1, NORBS110:                 DO K = 1, NORBS
115:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(I,J)*ROTATION(I,K)111:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(I,J)*ROTATION(I,K)
116:                 END DO112:                 END DO
117:             END DO113:             END DO
118:             !PRINT *,'ORBITAL',I,',ORBVAR=',ORBVAR(I),',DIPOLES=',ORBDIPOLES(1:3,I)114:             !PRINT *,'ORBITAL',I,',ORBVAR=',ORBVAR(I),',DIPOLES=',ORBDIPOLES(1:3,I)
 115: 
119:         END DO116:         END DO
120: 117: 
121:         !PRINT *,'ORBVARS=',ORBVAR118:         !PRINT *,'ORBVARS=',ORBVAR
122:         !PRINT *,'ORBVARSEXP=',ORBVAR**ORBVAREXPONENT119:         !PRINT *,'ORBVARSEXP=',ORBVAR**ORBVAREXPONENT
123: 120: 
124:         LOCALITY = SUM(ORBVAR ** ORBVAREXPONENT)121:         LOCALITY = SUM(ORBVAR ** ORBVAREXPONENT)
125: 122: 
126:         !PRINT *,'LOCALITY=',LOCALITY123:         !PRINT *,'LOCALITY=',LOCALITY
127: 124: 
128:         IF (GTEST) THEN125:         IF (GTEST) THEN
234:         DO TRIIND = 1, NROTS231:         DO TRIIND = 1, NROTS
235:             CALL DECODE_TRIIND(TRIIND, I, J)232:             CALL DECODE_TRIIND(TRIIND, I, J)
236:             IF (TRIIND.EQ.DERIVTRIIND) THEN233:             IF (TRIIND.EQ.DERIVTRIIND) THEN
237:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)234:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)
238:             ELSE235:             ELSE
239:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)236:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)
240:             END IF237:             END IF
241:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT)238:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT)
242:         END DO239:         END DO
243: 240: 
 241:         ! This appears to be required- investigate why. Something confusing in
 242:         ! fortran array indexing.
 243:         DERIV_ROTATION = TRANSPOSE(DERIV_ROTATION)
244:         !DO TRIIND = 1, NORBS244:         !DO TRIIND = 1, NORBS
245:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS)245:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS)
246:         !END DO246:         !END DO
247: 247: 
248:     END SUBROUTINE GET_ROTATION_DERIVATIVE248:     END SUBROUTINE GET_ROTATION_DERIVATIVE
249: 249: 
250:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT)250:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT)
251: 251: 
252:         ! Obtains Givens rotation between orbitals I and J with angle theta.252:         ! Obtains Givens rotation between orbitals I and J with angle theta.
253: 253: 
304:     END SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION304:     END SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION
305: 305: 
306:     SUBROUTINE GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)306:     SUBROUTINE GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)
307: 307: 
308:         ! Obtain the dipole moments of the orbital set created by the current rotation of the starting MOs.308:         ! Obtain the dipole moments of the orbital set created by the current rotation of the starting MOs.
309: 309: 
310:         USE COMMONS, ONLY: NORBS, DIPINTS310:         USE COMMONS, ONLY: NORBS, DIPINTS
311: 311: 
312:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS, NORBS)312:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS, NORBS)
313:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS)313:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS)
314:         ! DOUBLE PRECISION :: TEMPCURRENTDIPOLES(NORBS,NORBS) 
315: 314: 
316:         INTEGER :: CURRENT_ORB, ORIG_ORBA, ORIG_ORBB, DIRECTION315:         INTEGER :: CURRENT_ORB, ORIG_ORBA, ORIG_ORBB, DIRECTION
317: 316: 
318:         ORBDIPOLES(1:3,1:NORBS) = 0.0317:         ORBDIPOLES(1:3,1:NORBS) = 0.0
319: 318: 
320:         DO DIRECTION = 1,3319:         DO DIRECTION = 1,3
321:             DO CURRENT_ORB = 1, NORBS320:             DO CURRENT_ORB = 1, NORBS
322:                 DO ORIG_ORBA = 1, NORBS321:                 DO ORIG_ORBA = 1, NORBS
323:                     DO ORIG_ORBB = 1, NORBS322:                     DO ORIG_ORBB = 1, NORBS
324:                         ORBDIPOLES(DIRECTION,CURRENT_ORB) = ORBDIPOLES(DIRECTION,CURRENT_ORB) + &323:                         ORBDIPOLES(DIRECTION,CURRENT_ORB) = ORBDIPOLES(DIRECTION,CURRENT_ORB) + &
325:                             ROTATION(ORIG_ORBA,CURRENT_ORB) * ROTATION(ORIG_ORBB,CURRENT_ORB) * &324:                             ROTATION(CURRENT_ORB,ORIG_ORBA) * ROTATION(CURRENT_ORB,ORIG_ORBB) * &
326:                             DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB)325:                             DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB)
327:                     END DO326:                     END DO
328:                 END DO327:                 END DO
329:             END DO328:             END DO
330:         END DO329:         END DO
331: 330: 
332:         ! Comparison using matrix multiplication. 
333:         !DO DIRECTION = 1,3 
334:         !    PRINT *, 'Prev and matmul dipoles, direction', DIRECTION, '.' 
335:         !    TEMPCURRENTDIPOLES = MATMUL(MATMUL(TRANSPOSE(ROTATION),DIPINTS(DIRECTION,1:NORBS,1:NORBS)),ROTATION) 
336:         !    DO CURRENT_ORB = 1, NORBS 
337:         !        PRINT *, ORBDIPOLES(DIRECTION, CURRENT_ORB), TEMPCURRENTDIPOLES(CURRENT_ORB, CURRENT_ORB) 
338:         !    END DO 
339:         !END DO 
340:  
341:     END SUBROUTINE GET_CURRENT_ORB_DIPOLES331:     END SUBROUTINE GET_CURRENT_ORB_DIPOLES
342: 332: 
343:     SUBROUTINE DECODE_TRIIND(TRIIND, I, J)333:     SUBROUTINE DECODE_TRIIND(TRIIND, I, J)
344: 334: 
345:         ! Given a triangular index obtains the 2D indices corresponding to it.335:         ! Given a triangular index obtains the 2D indices corresponding to it.
346: 336: 
347:         USE COMMONS, ONLY: NORBS337:         USE COMMONS, ONLY: NORBS
348: 338: 
349:         INTEGER, INTENT(IN) :: TRIIND339:         INTEGER, INTENT(IN) :: TRIIND
350:         INTEGER, INTENT(OUT) :: I, J340:         INTEGER, INTENT(OUT) :: I, J
386:         DOUBLE PRECISION, INTENT(IN) :: CURRENT_ROTATION(NORBS, NORBS), CURRENT_DIP(3,NORBS)376:         DOUBLE PRECISION, INTENT(IN) :: CURRENT_ROTATION(NORBS, NORBS), CURRENT_DIP(3,NORBS)
387: 377: 
388:         DOUBLE PRECISION, INTENT(OUT) :: INDIVIDUAL_DERIVATIVES(NORBS,NORBS)378:         DOUBLE PRECISION, INTENT(OUT) :: INDIVIDUAL_DERIVATIVES(NORBS,NORBS)
389: 379: 
390:         INTEGER :: ORIGORB, NEWORB380:         INTEGER :: ORIGORB, NEWORB
391: 381: 
392:         !PRINT *,'CURRENTROTATION=',CURRENT_ROTATION382:         !PRINT *,'CURRENTROTATION=',CURRENT_ROTATION
393: 383: 
394:         DO NEWORB = 1, NORBS384:         DO NEWORB = 1, NORBS
395:             DO ORIGORB = 1, NORBS385:             DO ORIGORB = 1, NORBS
396:                 CALL GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEWORB, ORIGORB, INDIVIDUAL_DERIVATIVES(ORIGORB,NEWORB))386:                 CALL GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEWORB, ORIGORB, INDIVIDUAL_DERIVATIVES(NEWORB,ORIGORB))
397:             END DO387:             END DO
398:         END DO388:         END DO
399: 389: 
400:     END SUBROUTINE GET_INDIVIDUAL_DERIVATIVES390:     END SUBROUTINE GET_INDIVIDUAL_DERIVATIVES
401: 391: 
402:     SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEW_ORB, ORIG_ORB, DERIVATIVE)392:     SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEW_ORB, ORIG_ORB, DERIVATIVE)
403: 393: 
404:         ! Obtains the derivative of a single orbital at the current rotation's orbital394:         ! Obtains the derivative of a single orbital at the current rotation's orbital
405:         ! variance w.r.t. to the coefficient an original (unrotated) orbital.395:         ! variance w.r.t. to the coefficient an original (unrotated) orbital.
406: 396: 
412: 402: 
413:         INTEGER,INTENT(IN) :: NEW_ORB, ORIG_ORB403:         INTEGER,INTENT(IN) :: NEW_ORB, ORIG_ORB
414: 404: 
415:         INTEGER :: OTHER_ORIG_ORB405:         INTEGER :: OTHER_ORIG_ORB
416: 406: 
417:         DOUBLE PRECISION :: CONTRIB407:         DOUBLE PRECISION :: CONTRIB
418: 408: 
419:         DERIVATIVE = 0.0409:         DERIVATIVE = 0.0
420: 410: 
421:         DO OTHER_ORIG_ORB = 1, NORBS411:         DO OTHER_ORIG_ORB = 1, NORBS
422:             CONTRIB = 2 * CURRENT_ROTATION(OTHER_ORIG_ORB,NEW_ORB)*R2INTS(ORIG_ORB, OTHER_ORIG_ORB)412:             CONTRIB = 2 * CURRENT_ROTATION(NEW_ORB,OTHER_ORIG_ORB)*R2INTS(ORIG_ORB, OTHER_ORIG_ORB)
423: 413: 
424:             DERIVATIVE = DERIVATIVE + CONTRIB414:             DERIVATIVE = DERIVATIVE + CONTRIB
425: 415: 
426: 416: 
427:             CONTRIB = 4 * CURRENT_ROTATION(OTHER_ORIG_ORB,NEW_ORB)*&417:             CONTRIB = 4 * CURRENT_ROTATION(NEW_ORB,OTHER_ORIG_ORB)*&
428:                         (CURRENT_DIP(1,NEW_ORB)*DIPINTS(1,ORIG_ORB,OTHER_ORIG_ORB)&418:                         (CURRENT_DIP(1,NEW_ORB)*DIPINTS(1,ORIG_ORB,OTHER_ORIG_ORB)&
429:                         +CURRENT_DIP(2,NEW_ORB)*DIPINTS(2,ORIG_ORB,OTHER_ORIG_ORB)&419:                         +CURRENT_DIP(2,NEW_ORB)*DIPINTS(2,ORIG_ORB,OTHER_ORIG_ORB)&
430:                         +CURRENT_DIP(3,NEW_ORB)*DIPINTS(3,ORIG_ORB,OTHER_ORIG_ORB))420:                         +CURRENT_DIP(3,NEW_ORB)*DIPINTS(3,ORIG_ORB,OTHER_ORIG_ORB))
431: 421: 
432:             DERIVATIVE = DERIVATIVE - CONTRIB422:             DERIVATIVE = DERIVATIVE - CONTRIB
433: 423: 
434:         END DO424:         END DO
435: 425: 
436:     END SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE426:     END SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE
437: 427: 
444: 434: 
445:         DOUBLE PRECISION, INTENT(IN) :: IND_ORBVAR_DERIV(NORBS,NORBS), ORB_VAR(NORBS)435:         DOUBLE PRECISION, INTENT(IN) :: IND_ORBVAR_DERIV(NORBS,NORBS), ORB_VAR(NORBS)
446:         DOUBLE PRECISION, INTENT(OUT) :: PENALTY_DERIVATIVES(NORBS, NORBS)436:         DOUBLE PRECISION, INTENT(OUT) :: PENALTY_DERIVATIVES(NORBS, NORBS)
447: 437: 
448:         INTEGER :: ORIG_ORB, NEW_ORB438:         INTEGER :: ORIG_ORB, NEW_ORB
449: 439: 
450:         PENALTY_DERIVATIVES = 0.0440:         PENALTY_DERIVATIVES = 0.0
451: 441: 
452:         DO ORIG_ORB = 1, NORBS442:         DO ORIG_ORB = 1, NORBS
453:             DO NEW_ORB = 1, NORBS443:             DO NEW_ORB = 1, NORBS
454:                 PENALTY_DERIVATIVES(ORIG_ORB,NEW_ORB) = &444:                 PENALTY_DERIVATIVES(ORIG_ORB, NEW_ORB) = &
455:                     ORBVAREXPONENT * IND_ORBVAR_DERIV(ORIG_ORB,NEW_ORB) * ORB_VAR(NEW_ORB) ** (ORBVAREXPONENT-1)445:                     ORBVAREXPONENT * IND_ORBVAR_DERIV(NEW_ORB,ORIG_ORB) * ORB_VAR(NEW_ORB) ** (ORBVAREXPONENT-1)
456:             END DO446:             END DO
457:         END DO447:         END DO
458: 448: 
459:     END SUBROUTINE GET_PENALTY_DERIVATIVES449:     END SUBROUTINE GET_PENALTY_DERIVATIVES
460: 450: 
461:     SUBROUTINE ORBITALS_FINISH()451:     SUBROUTINE ORBITALS_FINISH()
462: 452: 
463:         ! Deallocate arrays in commons.453:         ! Deallocate arrays in commons.
464: 454: 
465:         USE COMMONS, ONLY: R2INTS, DIPINTS455:         USE COMMONS, ONLY: R2INTS, DIPINTS


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0