hdiff output

r32799/orbitals.f90 2017-06-15 15:30:14.621970010 +0100 r32798/orbitals.f90 2017-06-15 15:30:15.185977250 +0100
 82:         !PRINT '(A,I10,A)','NORBS=',NORBS 82:         !PRINT '(A,I10,A)','NORBS=',NORBS
 83:  83: 
 84:         DO I = 1, NORBS 84:         DO I = 1, NORBS
 85:            READ(INT_FILE_UNIT,*) INTS(I,1:NORBS)  85:            READ(INT_FILE_UNIT,*) INTS(I,1:NORBS) 
 86:            !PRINT *,INTS(I,1:NORBS) 86:            !PRINT *,INTS(I,1:NORBS)
 87:         END DO 87:         END DO
 88:  88: 
 89:         CLOSE(INT_FILE_UNIT) 89:         CLOSE(INT_FILE_UNIT)
 90:     END SUBROUTINE READ_INTEGRALS 90:     END SUBROUTINE READ_INTEGRALS
 91:  91: 
 92:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST, SECT) 92:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST)
 93:  93: 
 94:         ! Obtains the Power of the Orbital Variance locality measure for 94:         ! Obtains the Power of the Orbital Variance locality measure for
 95:         ! a given position. If GTEST is set also returns the gradient. 95:         ! a given position. If GTEST is set also returns the gradient.
 96:  96: 
 97:         ! Note that for intermediates dependent upon coefficient of current 97:         ! Note that for intermediates dependent upon coefficient of current
 98:         ! (rotated) orbital set in the original MO basis we use indexing 98:         ! (rotated) orbital set in the original MO basis we use indexing
 99:         ! of (ORIG_ORB,NEW_ORB) in all cases. 99:         ! of (ORIG_ORB,NEW_ORB) in all cases.
100: 100: 
101:         USE COMMONS, ONLY: NROTS, NORBS, R2INTS, DIPINTS, ORBVAREXPONENT101:         USE COMMONS, ONLY: NROTS, NORBS, R2INTS, DIPINTS, ORBVAREXPONENT
102: 102: 
103:         DOUBLE PRECISION, INTENT(INOUT) :: X(NROTS)103:         DOUBLE PRECISION, INTENT(IN) :: X(NROTS)
104:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(NROTS), LOCALITY104:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(NROTS), LOCALITY
105:         LOGICAL, INTENT(IN) :: GTEST,SECT105:         LOGICAL, INTENT(IN) :: GTEST
106: 106: 
107:         DOUBLE PRECISION :: ROTATION(NORBS,NORBS), ORBVAR(NORBS)107:         DOUBLE PRECISION :: ROTATION(NORBS,NORBS), ORBVAR(NORBS)
108:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS)108:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS)
109:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS)109:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS)
110:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS)110:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS)
111: 111: 
112:         INTEGER :: I,J,K112:         INTEGER :: I,J,K
113: 113: 
114:         CALL CHECK_COORDINATES(X) 
115:  
116:         CALL GET_ROTATION(X, ROTATION)114:         CALL GET_ROTATION(X, ROTATION)
117: 115: 
118:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)116:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)
119:         DO I = 1, NORBS117:         DO I = 1, NORBS
120:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2)118:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2)
121:             !PRINT*,ORBVAR(I)119:             !PRINT*,ORBVAR(I)
122:             DO J = 1, NORBS120:             DO J = 1, NORBS
123:                 DO K = 1, NORBS121:                 DO K = 1, NORBS
124:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(J,I)*ROTATION(K,I)122:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(J,I)*ROTATION(K,I)
125:                 END DO123:                 END DO
164:                 !PRINT *,'OVERALL_PROD='162:                 !PRINT *,'OVERALL_PROD='
165:                 !DO J = 1,NORBS163:                 !DO J = 1,NORBS
166:                 !    PRINT *, PROD_DERIVS(J,:)164:                 !    PRINT *, PROD_DERIVS(J,:)
167:                 !END DO165:                 !END DO
168:                 GRAD(I) = SUM(PROD_DERIVS)166:                 GRAD(I) = SUM(PROD_DERIVS)
169:                 !PRINT *,'GRAD',I,'=',GRAD(I)167:                 !PRINT *,'GRAD',I,'=',GRAD(I)
170:             END DO168:             END DO
171: 169: 
172:         END IF170:         END IF
173: 171: 
174:         IF (SECT) THEN 
175:             ! If haven't calculated gradient as well need to calculate additional intermediates. 
176:             IF (.NOT.GTEST) THEN 
177:                 CALL GET_INDIVIDUAL_DERIVATIVES(ROTATION, ORBDIPOLES, IND_ORBVAR_DERIV) 
178:                 CALL GET_PENALTY_DERIVATIVES(IND_ORBVAR_DERIV, ORBVAR, PENALTY_DERIV) 
179:             END IF 
180:  
181:             CALL CALC_HESSIAN(X,ROTATION,ORBVAR,ORBDIPOLES,IND_ORBVAR_DERIV,PENALTY_DERIV) 
182:  
183:         END IF 
184:         !PRINT*,"COORDS=",X,",GRAD=",GRAD 
185:     END SUBROUTINE GET_ORBITAL_LOCALITY172:     END SUBROUTINE GET_ORBITAL_LOCALITY
186: 173: 
187:     SUBROUTINE CHECK_COORDINATES(COORDS) 
188:  
189:         USE COMMONS, ONLY: NROTS 
190:  
191:         DOUBLE PRECISION, INTENT(INOUT) :: COORDS(1:NROTS) 
192:  
193:         INTEGER :: I, NDIFF 
194:  
195:         DOUBLE PRECISION, PARAMETER :: PI = 3.141592654D0 
196:         DOUBLE PRECISION :: COORDSSAVE(1:NROTS) 
197:  
198:         COORDSSAVE = COORDS 
199:  
200:         DO I=1,NROTS 
201:             IF (COORDS(I) > PI) THEN 
202:                 NDIFF = CEILING((COORDS(I)+PI)/(2*PI)) 
203:                 COORDS(I) = COORDS(I) - 2*PI*(NDIFF-1) 
204: !                PRINT*,COORDSSAVE(I),'->',COORDS(I),SIN(COORDSSAVE(I))-SIN(COORDS(I)) 
205:             ELSE IF (COORDS(I) < -PI) THEN 
206:                 NDIFF = CEILING((-COORDS(I)+PI)/(2*PI)) 
207:                 COORDS(I) = COORDS(I) + 2*PI*(NDIFF-1) 
208: !                PRINT*,COORDSSAVE(I),'->',COORDS(I),SIN(COORDSSAVE(I))-SIN(COORDS(I)) 
209:             END IF 
210:         END DO 
211:  
212:     END SUBROUTINE CHECK_COORDINATES 
213:  
214:     SUBROUTINE ELEMENTWISE_MULTIPLICATION(MATA,MATB,MATC)174:     SUBROUTINE ELEMENTWISE_MULTIPLICATION(MATA,MATB,MATC)
215: 175: 
216:         USE COMMONS, ONLY: NORBS176:         USE COMMONS, ONLY: NORBS
217: 177: 
218:         DOUBLE PRECISION, INTENT(IN) :: MATA(NORBS,NORBS), MATB(NORBS,NORBS)178:         DOUBLE PRECISION, INTENT(IN) :: MATA(NORBS,NORBS), MATB(NORBS,NORBS)
219: 179: 
220:         DOUBLE PRECISION, INTENT(OUT) :: MATC(NORBS, NORBS)180:         DOUBLE PRECISION, INTENT(OUT) :: MATC(NORBS, NORBS)
221: 181: 
222:         INTEGER :: I, J182:         INTEGER :: I, J
223: 183: 
280: 240: 
281:         DO TRIIND = 1, NROTS241:         DO TRIIND = 1, NROTS
282:             CALL DECODE_TRIIND(TRIIND, I, J)242:             CALL DECODE_TRIIND(TRIIND, I, J)
283:             IF (TRIIND.EQ.DERIVTRIIND) THEN243:             IF (TRIIND.EQ.DERIVTRIIND) THEN
284:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)244:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)
285:             ELSE245:             ELSE
286:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)246:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT)
287:             END IF247:             END IF
288:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT)248:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT)
289:         END DO249:         END DO
290:  
291:         !DO TRIIND = 1, NORBS250:         !DO TRIIND = 1, NORBS
292:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS)251:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS)
293:         !END DO252:         !END DO
294: 253: 
295:     END SUBROUTINE GET_ROTATION_DERIVATIVE254:     END SUBROUTINE GET_ROTATION_DERIVATIVE
296: 255: 
297:     SUBROUTINE GET_ROTATION_SECOND_DERIVATIVE(COORDS, DERIVTRIINDA, DERIVTRIINDB, SECDERIV_ROTATION) 
298:  
299:         ! Obtain the derivative of the rotation matrix with respect to a single 
300:         ! rotation angle theta at the current coordinates. 
301:  
302:         USE COMMONS, ONLY: NROTS, NORBS 
303:  
304:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
305:         INTEGER, INTENT(IN) :: DERIVTRIINDA, DERIVTRIINDB 
306:  
307:         DOUBLE PRECISION, INTENT(OUT) :: SECDERIV_ROTATION(NORBS,NORBS) 
308:  
309:         DOUBLE PRECISION :: NEXT_ROT(NORBS,NORBS) 
310:         INTEGER :: TRIIND, I, J 
311:  
312:         SECDERIV_ROTATION(1:NORBS,1:NORBS) = 0.0 
313:  
314:         DO TRIIND = 1, NORBS 
315:             SECDERIV_ROTATION(TRIIND, TRIIND) = 1.0 
316:         END DO 
317:  
318:         DO TRIIND = 1, NROTS 
319:             CALL DECODE_TRIIND(TRIIND, I, J) 
320:             IF (TRIIND.EQ.DERIVTRIINDA) THEN 
321:                 IF (DERIVTRIINDA.EQ.DERIVTRIINDB) THEN 
322:                     CALL GET_SECOND_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
323:                 ELSE 
324:                     CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
325:                 END IF 
326:             ELSE IF (TRIIND.EQ.DERIVTRIINDB) THEN 
327:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
328:             ELSE 
329:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
330:             END IF 
331:             SECDERIV_ROTATION = MATMUL(SECDERIV_ROTATION,NEXT_ROT) 
332:         END DO 
333:  
334:         !DO TRIIND = 1, NORBS 
335:         !    PRINT *,TRIIND,'DERIV_ROT=',SECDERIV_ROTATION(TRIIND,1:NORBS) 
336:         !END DO 
337:  
338:     END SUBROUTINE GET_ROTATION_SECOND_DERIVATIVE 
339:  
340:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT)256:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT)
341: 257: 
342:         ! Obtains Givens rotation between orbitals I and J with rotation angle theta.258:         ! Obtains Givens rotation between orbitals I and J with angle theta.
343: 259: 
344:         USE COMMONS, ONLY: NORBS260:         USE COMMONS, ONLY: NORBS
345: 261: 
346:         INTEGER, INTENT(IN) :: I, J262:         INTEGER, INTENT(IN) :: I, J
347:         DOUBLE PRECISION, INTENT(IN) :: THETA263:         DOUBLE PRECISION, INTENT(IN) :: THETA
348:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS)264:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS)
349: 265: 
350:         INTEGER :: VAL266:         INTEGER :: VAL
351: 267: 
352:         MAT(1:NORBS,1:NORBS) = 0.0268:         MAT(1:NORBS,1:NORBS) = 0.0
362: 278: 
363:         !PRINT *, 'GIVENSROT',I,'-',J,'='279:         !PRINT *, 'GIVENSROT',I,'-',J,'='
364:         !DO VAL=1,NORBS280:         !DO VAL=1,NORBS
365:         !    PRINT *,MAT(VAL,1:NORBS)281:         !    PRINT *,MAT(VAL,1:NORBS)
366:         !END DO282:         !END DO
367: 283: 
368:     END SUBROUTINE GET_GIVENS_ROTATION284:     END SUBROUTINE GET_GIVENS_ROTATION
369: 285: 
370:     SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION(I, J, THETA, MAT)286:     SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION(I, J, THETA, MAT)
371: 287: 
372:         ! Obtains the first derivative of the Givens rotation between orbitals I and J288:         ! Obtains the derivative of the Givens rotation between orbitals I and J with angle
373:         !  wrt the rotation angle at angle theta as a matrix.289:         ! theta as a matrix wrt the rotation angle theta.
374: 290: 
375:         USE COMMONS, ONLY: NORBS291:         USE COMMONS, ONLY: NORBS
376: 292: 
377:         INTEGER, INTENT(IN) :: I, J293:         INTEGER, INTENT(IN) :: I, J
378:         DOUBLE PRECISION, INTENT(IN) :: THETA294:         DOUBLE PRECISION, INTENT(IN) :: THETA
379:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS)295:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS)
380: 296: 
381:         INTEGER :: VAL297:         INTEGER :: VAL
382: 298: 
383:         MAT(1:NORBS,1:NORBS) = 0.0299:         MAT(1:NORBS,1:NORBS) = 0.0
386:         MAT(J,J) = -SIN(THETA)302:         MAT(J,J) = -SIN(THETA)
387:         MAT(I,J) = -COS(THETA)303:         MAT(I,J) = -COS(THETA)
388:         MAT(J,I) = COS(THETA)304:         MAT(J,I) = COS(THETA)
389:         !PRINT *, 'DERIVGIVENSROT',I,'-',J,'='305:         !PRINT *, 'DERIVGIVENSROT',I,'-',J,'='
390:         !DO VAL=1,NORBS306:         !DO VAL=1,NORBS
391:         !    PRINT *,MAT(VAL,1:NORBS)307:         !    PRINT *,MAT(VAL,1:NORBS)
392:         !END DO308:         !END DO
393: 309: 
394:     END SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION310:     END SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION
395: 311: 
396:     SUBROUTINE GET_SECOND_DERIVATIVE_GIVENS_ROTATION(I, J, THETA, MAT) 
397:  
398:         ! Obtains the second derivative of the Givens rotation between orbitals I and J 
399:         !  wrt the rotation angle at angle theta as a matrix. 
400:  
401:         USE COMMONS, ONLY: NORBS 
402:  
403:         INTEGER, INTENT(IN) :: I, J 
404:         DOUBLE PRECISION, INTENT(IN) :: THETA 
405:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS) 
406:  
407:         INTEGER :: VAL 
408:  
409:         MAT(1:NORBS,1:NORBS) = 0.0 
410:  
411:         MAT(I,I) = -COS(THETA) 
412:         MAT(J,J) = -COS(THETA) 
413:         MAT(I,J) = SIN(THETA) 
414:         MAT(J,I) = -SIN(THETA) 
415:         !PRINT *, 'SECONDDERIVGIVENSROT',I,'-',J,'=' 
416:         !DO VAL=1,NORBS 
417:         !    PRINT *,MAT(VAL,1:NORBS) 
418:         !END DO 
419:  
420:     END SUBROUTINE GET_SECOND_DERIVATIVE_GIVENS_ROTATION 
421:  
422:     SUBROUTINE GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)312:     SUBROUTINE GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES)
423: 313: 
424:         ! Obtain the dipole moments of the orbital set created by the current rotation of the starting MOs.314:         ! Obtain the dipole moments of the orbital set created by the current rotation of the starting MOs.
425: 315: 
426:         USE COMMONS, ONLY: NORBS, DIPINTS316:         USE COMMONS, ONLY: NORBS, DIPINTS
427: 317: 
428:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS, NORBS)318:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS, NORBS)
429:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS)319:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS)
430:         ! DOUBLE PRECISION :: TEMPCURRENTDIPOLES(NORBS,NORBS)320:         ! DOUBLE PRECISION :: TEMPCURRENTDIPOLES(NORBS,NORBS)
431: 321: 
567: 457: 
568:         DO ORIG_ORB = 1, NORBS458:         DO ORIG_ORB = 1, NORBS
569:             DO NEW_ORB = 1, NORBS459:             DO NEW_ORB = 1, NORBS
570:                 PENALTY_DERIVATIVES(ORIG_ORB,NEW_ORB) = &460:                 PENALTY_DERIVATIVES(ORIG_ORB,NEW_ORB) = &
571:                     ORBVAREXPONENT * IND_ORBVAR_DERIV(ORIG_ORB,NEW_ORB) * ORB_VAR(NEW_ORB) ** (ORBVAREXPONENT-1)461:                     ORBVAREXPONENT * IND_ORBVAR_DERIV(ORIG_ORB,NEW_ORB) * ORB_VAR(NEW_ORB) ** (ORBVAREXPONENT-1)
572:             END DO462:             END DO
573:         END DO463:         END DO
574: 464: 
575:     END SUBROUTINE GET_PENALTY_DERIVATIVES465:     END SUBROUTINE GET_PENALTY_DERIVATIVES
576: 466: 
577:     SUBROUTINE CALC_HESSIAN(COORDS,ROTATION,ORBVAR,ORBDIPOLES,IND_ORBVAR_DERIV,PENALTY_DERIV) 
578:  
579:         USE COMMONS, ONLY: NROTS, NORBS 
580:         USE MODHESS 
581:  
582:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
583:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS,NORBS), ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 
584:         DOUBLE PRECISION, INTENT(IN) :: ORBVAR(NORBS), PENALTY_DERIV(NORBS,NORBS) 
585:  
586:         DOUBLE PRECISION :: IND_ORBVAR_SECOND_DERIV(NORBS,NORBS,NORBS) 
587:         DOUBLE PRECISION :: PENALTY_SECOND_DERIV(NORBS,NORBS,NORBS), TEMP 
588:         DOUBLE PRECISION :: ROTMATA(NORBS,NORBS), ROTMATB(NORBS,NORBS) 
589:  
590:         INTEGER :: TRIINDA, TRIINDB, NEW_ORB, ORIG_ORBA, ORIG_ORBB 
591:  
592:         !IF (.NOT.ALLOCATED(HESS)) THEN 
593:         !    PRINT *, 'Allocating Hessian with size ',NROTS,'x',NROTS 
594:         !    ALLOCATE(HESS(NROTS,NROTS)) 
595:         !END IF 
596:         !PRINT*,'hessian bounds:',LBOUND(HESS,DIM=1),UBOUND(HESS,DIM=1),LBOUND(HESS,DIM=2),UBOUND(HESS,DIM=2) 
597:         !PRINT*,'hessian size:',SIZE(HESS),SHAPE(HESS) 
598:         ! Gradually work our way up the chain rule... 
599:         CALL CALC_ORBVAR_SECOND_DERIV(ROTATION,ORBDIPOLES,IND_ORBVAR_SECOND_DERIV) 
600:  
601:         CALL CALC_PENALTY_SECOND_DERIV(ORBVAR,IND_ORBVAR_DERIV,IND_ORBVAR_SECOND_DERIV,PENALTY_SECOND_DERIV) 
602:  
603:  
604:         !PRINT*,NROTS,NORBS 
605:  
606:         !PRINT*,'IND_ORBVAR_SECOND_DERIV=' 
607:         !DO TRIINDA=1,NORBS 
608:         !    DO TRIINDB=1,NORBS 
609:         !        PRINT*,TRIINDA,TRIINDB,IND_ORBVAR_SECOND_DERIV(TRIINDA,TRIINDB,1:NORBS) 
610:         !    END DO 
611:         !END DO 
612:         !PRINT*,'PENALTY_SECOND_DERIV=' 
613:         !DO TRIINDA=1,NORBS 
614:         !    DO TRIINDB=1,NORBS 
615:         !        PRINT*,TRIINDA,TRIINDB,PENALTY_SECOND_DERIV(TRIINDA,TRIINDB,1:NORBS) 
616:         !    END DO 
617:         !END DO 
618:  
619:         HESS(1:NROTS,1:NROTS) = 0.0 
620:  
621:         ! Now calculate actual hessian! 
622:         DO TRIINDA=1,NROTS 
623:             DO TRIINDB=1,NROTS 
624:                 ! First deal with term proportional to penalty first deriv wrt coeffs 
625:                 CALL GET_ROTATION_SECOND_DERIVATIVE(COORDS,TRIINDA,TRIINDB,ROTMATA) 
626:  
627:                 ! Use ROTMATB as scratch space. 
628:                 CALL ELEMENTWISE_MULTIPLICATION(ROTMATA(1:NORBS,1:NORBS),PENALTY_DERIV(1:NORBS,1:NORBS),ROTMATB) 
629:                 !PRINT*,'PRODMAT SECDERIV=' 
630:                 !DO NEW_ORB = 1,NORBS 
631:                 !    PRINT*,ROTMATB(NEW_ORB,1:NORBS) 
632:                 !END DO 
633:  
634:                 HESS(TRIINDA,TRIINDB) = SUM(ROTMATB) 
635:                 ! Now deal with term proportional to penalty second deriv wrt coeffs. 
636:                 CALL GET_ROTATION_DERIVATIVE(COORDS,TRIINDA,ROTMATA) 
637:                 CALL GET_ROTATION_DERIVATIVE(COORDS,TRIINDB,ROTMATB) 
638:                 DO NEW_ORB = 1,NORBS 
639:                     DO ORIG_ORBA = 1,NORBS 
640:                         DO ORIG_ORBB = 1,NORBS 
641:                             HESS(TRIINDA,TRIINDB) = HESS(TRIINDA,TRIINDB) + PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB)*ROTMATA(ORIG_ORBA,NEW_ORB)*ROTMATB(ORIG_ORBB,NEW_ORB) 
642:                         END DO 
643:                     END DO 
644:                 END DO 
645:             END DO 
646:         END DO 
647:  
648:     END SUBROUTINE CALC_HESSIAN 
649:  
650:     SUBROUTINE CALC_ORBVAR_SECOND_DERIV(ROTATION,ORBDIPOLES,IND_ORBVAR_SECOND_DERIV) 
651:  
652:         USE COMMONS, ONLY: NORBS, R2INTS, DIPINTS 
653:  
654:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS,NORBS), ORBDIPOLES(3,NORBS) 
655:         DOUBLE PRECISION, INTENT(OUT) :: IND_ORBVAR_SECOND_DERIV(NORBS,NORBS,NORBS) 
656:  
657:         DOUBLE PRECISION :: INTERMEDIATE_SUM(NORBS) 
658:  
659:         INTEGER :: ORIG_ORBA, ORIG_ORBB, NEW_ORB, DIRECTION 
660:  
661:         ! First set values according to [x2] term. 
662:         DO ORIG_ORBA = 1,NORBS 
663:             DO ORIG_ORBB = 1,NORBS 
664:                 IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,1:NORBS) = 2 * R2INTS(ORIG_ORBA,ORIG_ORBB) 
665:             END DO 
666:         END DO 
667:  
668:         DO NEW_ORB = 1, NORBS 
669:             DO DIRECTION = 1,3 
670:                 ! Calc \sum_i C_ip [x]_ij as for each j value (and given p). 
671:                 INTERMEDIATE_SUM(1:NORBS) = 0.0 
672:                 DO ORIG_ORBA = 1,NORBS 
673:                     DO ORIG_ORBB = 1,NORBS 
674:                         INTERMEDIATE_SUM(ORIG_ORBA) = INTERMEDIATE_SUM(ORIG_ORBA) + ROTATION(ORIG_ORBB,NEW_ORB) * DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB) 
675:                     END DO 
676:                 END DO 
677:                 !PRINT*,'DIRECTION=',DIRECTION 
678:                 !PRINT*,'DIPINTS=',DIPINTS(DIRECTION,:,:) 
679:                 !PRINT*,'INTERMEDIATE_SUM=',INTERMEDIATE_SUM 
680:                 DO ORIG_ORBA = 1,NORBS 
681:                     DO ORIG_ORBB = 1,NORBS 
682:                         IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) = IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) & 
683:                                             - 8*INTERMEDIATE_SUM(ORIG_ORBA)*INTERMEDIATE_SUM(ORIG_ORBB) & 
684:                                             - 4*ORBDIPOLES(DIRECTION,NEW_ORB)*DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB) 
685:                     END DO 
686:                 END DO 
687:  
688:             END DO 
689:         END DO 
690:  
691:     END SUBROUTINE CALC_ORBVAR_SECOND_DERIV 
692:  
693:     SUBROUTINE CALC_PENALTY_SECOND_DERIV(ORBVAR,IND_ORBVAR_DERIV,IND_ORBVAR_SECOND_DERIV,PENALTY_SECOND_DERIV) 
694:  
695:         USE COMMONS, ONLY: ORBVAREXPONENT, NORBS 
696:  
697:         DOUBLE PRECISION, INTENT(IN) :: ORBVAR(NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 
698:  
699:         DOUBLE PRECISION, INTENT(IN) :: IND_ORBVAR_SECOND_DERIV(NORBS,NORBS,NORBS) 
700:         DOUBLE PRECISION, INTENT(OUT) :: PENALTY_SECOND_DERIV(NORBS,NORBS,NORBS) 
701:  
702:         INTEGER :: NEW_ORB, ORIG_ORBA, ORIG_ORBB 
703:  
704:         DO NEW_ORB = 1,NORBS 
705:             DO ORIG_ORBA = 1,NORBS 
706:                 DO ORIG_ORBB = 1,NORBS 
707:                     PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) = & 
708:                                 ORBVAREXPONENT * IND_ORBVAR_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB)*(ORBVAR(NEW_ORB))**(ORBVAREXPONENT-1) 
709:                     IF (ORBVAREXPONENT.GT.1) THEN 
710:                         PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB) = PENALTY_SECOND_DERIV(ORIG_ORBA,ORIG_ORBB,NEW_ORB)& 
711:                                 +ORBVAREXPONENT*(ORBVAREXPONENT-1)*IND_ORBVAR_DERIV(ORIG_ORBA,NEW_ORB)*IND_ORBVAR_DERIV(ORIG_ORBB,NEW_ORB)*(ORBVAR(NEW_ORB))**(ORBVAREXPONENT-2) 
712:                     END IF 
713:                 END DO 
714:             END DO 
715:         END DO 
716:  
717:     END SUBROUTINE CALC_PENALTY_SECOND_DERIV 
718:  
719:     SUBROUTINE ORBITALS_FINISH()467:     SUBROUTINE ORBITALS_FINISH()
720: 468: 
721:         ! Deallocate arrays in commons.469:         ! Deallocate arrays in commons.
722: 470: 
723:         USE COMMONS, ONLY: R2INTS, DIPINTS471:         USE COMMONS, ONLY: R2INTS, DIPINTS
724: 472: 
725:         IF (ALLOCATED(R2INTS)) DEALLOCATE(R2INTS)473:         IF (ALLOCATED(R2INTS)) DEALLOCATE(R2INTS)
726:         IF (ALLOCATED(DIPINTS)) DEALLOCATE(DIPINTS)474:         IF (ALLOCATED(DIPINTS)) DEALLOCATE(DIPINTS)
727: 475: 
728:     END SUBROUTINE ORBITALS_FINISH476:     END SUBROUTINE ORBITALS_FINISH


r32799/potential.f90 2017-06-15 15:30:14.905973649 +0100 r32798/potential.f90 2017-06-15 15:30:15.465980845 +0100
425: ! !        IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN425: ! !        IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN
426: !             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)426: !             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
427: ! !        ENDIF427: ! !        ENDIF
428: !       ENDDO428: !       ENDDO
429: 429: 
430:    ELSE IF (MLP3T) THEN430:    ELSE IF (MLP3T) THEN
431:       CALL MLP3(X, GRAD, EREAL, GRADT, SECT)431:       CALL MLP3(X, GRAD, EREAL, GRADT, SECT)
432:    ELSEIF (MLQT) THEN432:    ELSEIF (MLQT) THEN
433:       CALL MLQ(X, GRAD, EREAL, GRADT, SECT)433:       CALL MLQ(X, GRAD, EREAL, GRADT, SECT)
434:    ELSE IF (ORBITALS) THEN434:    ELSE IF (ORBITALS) THEN
435:        CALL GET_ORBITAL_LOCALITY(X, GRAD, EREAL, GRADT,SECT)435:       CALL GET_ORBITAL_LOCALITY(X, GRAD, EREAL, GRADT)
436: !       DIFF=1.0D-4436: !       DIFF=1.0D-4
437: !       WRITE(MYUNIT, *) 'analytic and numerical gradients:'437: !       WRITE(MYUNIT, *) 'analytic and numerical gradients:'
438: !       DO J1=1, NATOMS438: !       DO J1=1, NATOMS
439: !          X(J1)=X(J1)+DIFF439: !          X(J1)=X(J1)+DIFF
440: !          CALL GET_ORBITAL_LOCALITY(X, GPLUS, EPLUS,.FALSE.,.FALSE.)440: !          CALL GET_ORBITAL_LOCALITY(X, GPLUS, EPLUS,.FALSE.)
441: !          X(J1)=X(J1)-2.0D0*DIFF441: !          X(J1)=X(J1)-2.0D0*DIFF
442: !          CALL GET_ORBITAL_LOCALITY(X, GMINUS, EMINUS,.FALSE.,.FALSE.)442: !          CALL GET_ORBITAL_LOCALITY(X, GMINUS, EMINUS,.FALSE.)
443: !          X(J1)=X(J1)+DIFF443: !          X(J1)=X(J1)+DIFF
444: !          IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN444: !          IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN
445: !             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)445: !             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
446: !          ENDIF446: !          ENDIF
447: !       ENDDO447: !       ENDDO
448: !       IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(NROTS,NROTS)) 
449: !       IF (.NOT.SECT) CALL GET_ORBITAL_LOCALITY(X, GRAD, EREAL, GRADT,.TRUE.) 
450: !       WRITE(MYUNIT,*) 'analytic and numerical second derivatives:' 
451: !       DO J1=1,NROTS 
452: !            X(J1)=X(J1)+DIFF 
453: !            CALL GET_ORBITAL_LOCALITY(X,GPLUS,EPLUS,.TRUE.,.FALSE.) 
454: !            X(J1)=X(J1)-2.0D0*DIFF 
455: !            CALL GET_ORBITAL_LOCALITY(X,GMINUS,EMINUS,.TRUE.,.FALSE.) 
456: !            X(J1)=X(J1)+DIFF 
457: !            DO J2=1,NROTS 
458: !               IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND. & 
459: !                    (ABS(100.0D0*(HESS(J1,J2)-(GPLUS(J2)-GMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN 
460: !                  WRITE(MYUNIT,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(GPLUS(J2)-GMINUS(J2))/(2.0D0*DIFF),'   X' 
461: !               ELSE 
462: !                  WRITE(MYUNIT,'(2I5,2F20.10,A)') J1,J2,HESS(J1,J2),(GPLUS(J2)-GMINUS(J2))/(2.0D0*DIFF) 
463: !               ENDIF 
464: !            ENDDO 
465: !       ENDDO 
466: !       STOP448: !       STOP
467:    ELSE IF (LJATT) THEN449:    ELSE IF (LJATT) THEN
468:       CALL LJ(X, GRAD, EREAL, GRADT, SECT)450:       CALL LJ(X, GRAD, EREAL, GRADT, SECT)
469:       CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)451:       CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)
470: 452: 
471:    ELSE IF (DFTBCT) THEN453:    ELSE IF (DFTBCT) THEN
472:       IF (.NOT.PERCOLATET) CALL RAD(X, GRAD, EREAL, GRADT)454:       IF (.NOT.PERCOLATET) CALL RAD(X, GRAD, EREAL, GRADT)
473:       CALL DFTBC(NATOMS, X, GRAD, EREAL, GRADT)455:       CALL DFTBC(NATOMS, X, GRAD, EREAL, GRADT)
474:       IF (FTEST) THEN456:       IF (FTEST) THEN
475:          RETURN457:          RETURN


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0