hdiff output

r31729/potential.f 2017-01-21 10:35:18.994250000 +0000 r31728/potential.f 2017-01-21 10:35:19.226250000 +0000
311:                   ! PRINT '(6G20.10)',HESS(1:NOPT,1:NOPT)311:                   ! PRINT '(6G20.10)',HESS(1:NOPT,1:NOPT)
312:                ENDIF312:                ENDIF
313:             ELSE313:             ELSE
314:                NZERO=0314:                NZERO=0
315:                ENERGY=0.0D0315:                ENERGY=0.0D0
316:                IF (GTEST) VNEW(1:NOPT)=0.0D0316:                IF (GTEST) VNEW(1:NOPT)=0.0D0
317:                IF (STEST) HESS(1:NOPT,1:NOPT)=0.0D0317:                IF (STEST) HESS(1:NOPT,1:NOPT)=0.0D0
318:                IF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'NIMET') THEN ! nickel(100) metal + hydrogen318:                IF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'NIMET') THEN ! nickel(100) metal + hydrogen
319:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units habr =1319:                   DUMMY2=1.0D0/(RPBETA/RPIMAGES)**2 ! atomic units habr =1
320:                   DO J1=1,RPIMAGES320:                   DO J1=1,RPIMAGES
321:                      CALL NIMETP(COORDS(RPDOF*(J1-1)+1:RPDOF*J1), VNEW(RPDOF*(J1-1)+1:RPDOF*J1), DUMMY1, GTEST, SSTEST)321:                      CALL NIMETP(COORDS(RPDOF*(J1-1)+1:RPDOF*J1), VNEW(RPDOF*(J1-1)+1:RPDOF*J1), DUMMY1, GTEST, STEST)
322:                      print*, ENERGY, 'EneTest'322:                      print*, ENERGY, 'EneTest'
323:                      ! ENERGY323:                      ! ENERGY
324:                      ! IF (GTEST.OR.STEST)324:                      ! IF (GTEST.OR.STEST)
325:                   ENDDO325:                   ENDDO
326: 326: 
327:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'SD') THEN ! Stillinger-David flexible water potential327:                ELSEIF (TRIM(ADJUSTL(RPSYSTEM)).EQ.'SD') THEN ! Stillinger-David flexible water potential
328:                   DUMMY2=1.0D0/(RPBETA*15.1787D0/RPIMAGES)**2 ! hbar = 15.1787 fs kcal / mol328:                   DUMMY2=1.0D0/(RPBETA*15.1787D0/RPIMAGES)**2 ! hbar = 15.1787 fs kcal / mol
329:                   DO J1=1,RPIMAGES329:                   DO J1=1,RPIMAGES
330:                      ENERGY=ENERGY+SDPOTENTIAL(COORDS(RPDOF*(J1-1)+1:RPDOF*J1))330:                      ENERGY=ENERGY+SDPOTENTIAL(COORDS(RPDOF*(J1-1)+1:RPDOF*J1))
331:                      IF (GTEST.OR.STEST) VNEW(RPDOF*(J1-1)+1:RPDOF*J1)=SDGRAD(COORDS(RPDOF*(J1-1)+1:RPDOF*J1))331:                      IF (GTEST.OR.STEST) VNEW(RPDOF*(J1-1)+1:RPDOF*J1)=SDGRAD(COORDS(RPDOF*(J1-1)+1:RPDOF*J1))
478:             478:             
479:          ELSE IF(MULTIPOTT) THEN479:          ELSE IF(MULTIPOTT) THEN
480:             CALL MULTIPOT_CALL(COORDS,VNEW,ENERGY,GTEST,SSTEST)480:             CALL MULTIPOT_CALL(COORDS,VNEW,ENERGY,GTEST,SSTEST)
481:             IF (PTEST) THEN481:             IF (PTEST) THEN
482:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '482:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
483:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '483:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
484:             ENDIF484:             ENDIF
485: 485: 
486:          ELSE IF(EX1DT) THEN486:          ELSE IF(EX1DT) THEN
487:             487:             
488:             CALL example1Dpotential(COORDS,VNEW,ENERGY,GTEST,SSTEST)488:             CALL example1Dpotential(COORDS,VNEW,ENERGY,GTEST,STEST)
489:             489:             
490:             IF (PTEST) THEN490:             IF (PTEST) THEN
491:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '491:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
492:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '492:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
493:             ENDIF493:             ENDIF
494: 494: 
495: 495: 
496:          ELSE IF (VARIABLES) THEN496:          ELSE IF (VARIABLES) THEN
497: 497: 
498:             IF (MLPVB3T) THEN498:             IF (MLPVB3T) THEN
499:                CALL MLPVB3(COORDS, VNEW, ENERGY, GTEST, SSTEST)499:                CALL MLPVB3(COORDS, VNEW, ENERGY, GTEST, STEST)
500: 500: 
501: !               DIFF=1.0D-4501: !               DIFF=1.0D-4
502: !               PRINT*,'analytic and numerical gradients:'502: !               PRINT*,'analytic and numerical gradients:'
503: !               IF (.NOT.(ALLOCATED(HESS))) ALLOCATE(HESS(NMLP,NMLP))503: !               IF (.NOT.(ALLOCATED(HESS))) ALLOCATE(HESS(NMLP,NMLP))
504: !               CALL MLPVB3(COORDS, VNEW, ENERGY, .TRUE., .TRUE.)504: !               CALL MLPVB3(COORDS, VNEW, ENERGY, .TRUE., .TRUE.)
505: !               PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS)505: !               PRINT '(A,I8)','SIZE(HESS)=',SIZE(HESS)
506: !               HESSDUM(1:NMLP,1:NMLP)=HESS(1:NMLP,1:NMLP)506: !               HESSDUM(1:NMLP,1:NMLP)=HESS(1:NMLP,1:NMLP)
507: !               DO J1=1,NATOMS507: !               DO J1=1,NATOMS
508: !                  COORDS(J1)=COORDS(J1)+DIFF508: !                  COORDS(J1)=COORDS(J1)+DIFF
509: !                  CALL MLPVB3(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)509: !                  CALL MLPVB3(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
528: !                     IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND. 528: !                     IF ((ABS(HESS(J1,J2)).NE.0.0D0).AND. 
529: !     &                  (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN529: !     &                  (ABS(100.0D0*(HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D0)) THEN
530: !                     WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'530: !                     WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),'   X'
531: !                     ELSE531: !                     ELSE
532: !                        WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)532: !                        WRITE(*,'(2I5,2G20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF)
533: !                     ENDIF533: !                     ENDIF
534: !                  ENDDO534: !                  ENDDO
535: !                ENDDO535: !                ENDDO
536: !                STOP536: !                STOP
537:             ELSEIF (MLPB3T) THEN537:             ELSEIF (MLPB3T) THEN
538:                CALL MLPB3(COORDS, VNEW, ENERGY, GTEST, SSTEST)538:                CALL MLPB3(COORDS, VNEW, ENERGY, GTEST, STEST)
539:             ELSEIF (MLP3T) THEN539:             ELSEIF (MLP3T) THEN
540:                CALL MLP3(COORDS, VNEW, ENERGY, GTEST, SSTEST)540:                CALL MLP3(COORDS, VNEW, ENERGY, GTEST, STEST)
541:             ELSEIF (MLQT) THEN541:             ELSEIF (MLQT) THEN
542:                CALL MLQ(COORDS, VNEW, ENERGY, GTEST, SSTEST)542:                CALL MLQ(COORDS, VNEW, ENERGY, GTEST, STEST)
543:             ELSE543:             ELSE
544:                CALL FUNCTIONAL( COORDS, VNEW, ENERGY, GTEST, SSTEST)544:                CALL FUNCTIONAL( COORDS, VNEW, ENERGY, GTEST, STEST)
545:             ENDIF545:             ENDIF
546:             IF (PTEST.OR.MLPPROB.OR.MLQPROB) THEN546:             IF (PTEST.OR.MLPPROB.OR.MLQPROB) THEN
547:                 WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '547:                 WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
548:                 WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '548:                 WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
549:             ENDIF549:             ENDIF
550: !           IF (MLPPROB) STOP550: !           IF (MLPPROB) STOP
551:             ! CALL CTEST(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST)551:             ! CALL CTEST(NATOMS, COORDS, VNEW, ENERGY, GTEST, STEST)
552:             ! CALL TWODFUNC(COORDS,VNEW,ENERGY,GTEST,SSTEST)552:             ! CALL TWODFUNC(COORDS,VNEW,ENERGY,GTEST,STEST)
553:             ! CALL MB(COORDS,VNEW,ENERGY,GTEST,SSTEST)553:             ! CALL MB(COORDS,VNEW,ENERGY,GTEST,STEST)
554:             ! CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,GTEST,SSTEST)554:             ! CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,GTEST,STEST)
555:             ! CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,GTEST,SSTEST)555:             ! CALL P4DIFF(NATOMS,COORDS,VNEW,ENERGY,PARAM1,PARAM2,GTEST,STEST)
556: !           CALL WATERMETHANE(COORDS,ENERGY)556: !           CALL WATERMETHANE(COORDS,ENERGY)
557: !           DIFF=1.0D-4557: !           DIFF=1.0D-4
558: !           DO J1=1,6558: !           DO J1=1,6
559: !              COORDS(J1)=COORDS(J1)+DIFF559: !              COORDS(J1)=COORDS(J1)+DIFF
560: !              CALL WATERMETHANE(COORDS,EPLUS)560: !              CALL WATERMETHANE(COORDS,EPLUS)
561: !              COORDS(J1)=COORDS(J1)-2.0D0*DIFF561: !              COORDS(J1)=COORDS(J1)-2.0D0*DIFF
562: !              CALL WATERMETHANE(COORDS,EMINUS)562: !              CALL WATERMETHANE(COORDS,EMINUS)
563: !              COORDS(J1)=COORDS(J1)+DIFF563: !              COORDS(J1)=COORDS(J1)+DIFF
564: !              VNEW(J1)=(EPLUS-EMINUS)/(2*DIFF)564: !              VNEW(J1)=(EPLUS-EMINUS)/(2*DIFF)
565: !           ENDDO565: !           ENDDO
567: 567: 
568:             IF (PHI4MODT) THEN568:             IF (PHI4MODT) THEN
569:                 CALL PHI4MODEL(COORDS,VNEW,ENERGY,GTEST,SSTEST)569:                 CALL PHI4MODEL(COORDS,VNEW,ENERGY,GTEST,SSTEST)
570:                 IF (PTEST) THEN570:                 IF (PTEST) THEN
571:                     WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '571:                     WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
572:                     WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '572:                     WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
573:                 ENDIF573:                 ENDIF
574: 574: 
575:          ELSE   IF (ONEDAPBCT) THEN575:          ELSE   IF (ONEDAPBCT) THEN
576: 576: 
577:                CALL ENERGY_1D_APBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)577:                CALL ENERGY_1D_APBC(COORDS,VNEW,ENERGY,GTEST,STEST)
578:                ! Debug tools578:                ! Debug tools
579:                ! 579:                ! 
580:                ! DIFF=1.0D-2580:                ! DIFF=1.0D-2
581:                ! PRINT*,'analytic and numerical gradients: NATOMS=',NATOMS581:                ! PRINT*,'analytic and numerical gradients: NATOMS=',NATOMS
582:                ! DO J1=1,NATOMS582:                ! DO J1=1,NATOMS
583:                ! COORDS(J1)=COORDS(J1)+DIFF583:                ! COORDS(J1)=COORDS(J1)+DIFF
584:                ! CALL ENERGY_1D_APBC(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)584:                ! CALL ENERGY_1D_APBC(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
585:                ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF585:                ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF
586:                ! CALL ENERGY_1D_APBC(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)586:                ! CALL ENERGY_1D_APBC(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)
587:                ! COORDS(J1)=COORDS(J1)+DIFF587:                ! COORDS(J1)=COORDS(J1)+DIFF
605:                ! ENDDO605:                ! ENDDO
606:                ! ENDDO606:                ! ENDDO
607: 607: 
608:                IF (PTEST) THEN608:                IF (PTEST) THEN
609:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '609:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
610:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '610:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
611:                ENDIF611:                ENDIF
612: 612: 
613:             ELSE IF(ONEDPBCT) THEN613:             ELSE IF(ONEDPBCT) THEN
614: 614: 
615:                CALL ENERGY_1D_PBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)615:                CALL ENERGY_1D_PBC(COORDS,VNEW,ENERGY,GTEST,STEST)
616: 616: 
617:                IF (PTEST) THEN617:                IF (PTEST) THEN
618:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '618:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
619:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '619:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
620:                ENDIF620:                ENDIF
621: 621: 
622: 622: 
623:             ELSE IF(INVTONEDPBCT) THEN623:             ELSE IF(INVTONEDPBCT) THEN
624: 624: 
625: !              CALL ENERGY_1D_PBC_INVT(COORDS,VNEW,ENERGY,GTEST,SSTEST)625: !              CALL ENERGY_1D_PBC_INVT(COORDS,VNEW,ENERGY,GTEST,STEST)
626: 626: 
627:                IF (PTEST) THEN627:                IF (PTEST) THEN
628:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '628:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
629:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '629:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
630:                ENDIF630:                ENDIF
631: 631: 
632:             ELSE IF(INVTTWODPBCT) THEN632:             ELSE IF(INVTTWODPBCT) THEN
633: 633: 
634:                CALL ENERGY_2D_PBC_INVT(COORDS,VNEW,ENERGY,GTEST,SSTEST)634:                CALL ENERGY_2D_PBC_INVT(COORDS,VNEW,ENERGY,GTEST,STEST)
635: 635: 
636:                IF (PTEST) THEN636:                IF (PTEST) THEN
637:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '637:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
638:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '638:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
639:                ENDIF639:                ENDIF
640: 640: 
641: 641: 
642:             ELSE IF(TWODAPBCT) THEN642:             ELSE IF(TWODAPBCT) THEN
643: 643: 
644:                IF(NATOMS.NE.(NONEDAPBC**2)) THEN644:                IF(NATOMS.NE.(NONEDAPBC**2)) THEN
645:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'645:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'
646:                ENDIF646:                ENDIF
647: 647: 
648:                CALL ENERGY_2D_APBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)648:                CALL ENERGY_2D_APBC(COORDS,VNEW,ENERGY,GTEST,STEST)
649: 649: 
650:                ! Debug tools650:                ! Debug tools
651:                ! DIFF=1.0D-2651:                ! DIFF=1.0D-2
652:                ! PRINT*,'analytic and numerical gradients: N_lattice, NATOMS=',(NONEDAPBC), NATOMS652:                ! PRINT*,'analytic and numerical gradients: N_lattice, NATOMS=',(NONEDAPBC), NATOMS
653:                ! DO J1=1,NATOMS653:                ! DO J1=1,NATOMS
654:                ! COORDS(J1)=COORDS(J1)+DIFF654:                ! COORDS(J1)=COORDS(J1)+DIFF
655:                ! CALL ENERGY_2D_APBC(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)655:                ! CALL ENERGY_2D_APBC(COORDS,VPLUS,EPLUS,.FALSE.,.FALSE.)
656:                ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF656:                ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF
657:                ! CALL ENERGY_2D_APBC(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)657:                ! CALL ENERGY_2D_APBC(COORDS,VMINUS,EMINUS,.FALSE.,.FALSE.)
658:                ! COORDS(J1)=COORDS(J1)+DIFF658:                ! COORDS(J1)=COORDS(J1)+DIFF
680:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '680:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
681:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '681:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
682:                ENDIF682:                ENDIF
683: 683: 
684:             ELSE IF(TWODPBCT) THEN684:             ELSE IF(TWODPBCT) THEN
685: 685: 
686:                IF(NATOMS.NE.(NONEDAPBC**2)) THEN686:                IF(NATOMS.NE.(NONEDAPBC**2)) THEN
687:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'687:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'
688:                ENDIF688:                ENDIF
689: 689: 
690:                CALL ENERGY_2D_PBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)690:                CALL ENERGY_2D_PBC(COORDS,VNEW,ENERGY,GTEST,STEST)
691: 691: 
692:                IF (PTEST) THEN692:                IF (PTEST) THEN
693:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '693:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
694:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '694:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
695:                ENDIF695:                ENDIF
696: 696: 
697: 697: 
698: 698: 
699: 699: 
700:             ELSE IF(THREEDAPBCT) THEN700:             ELSE IF(THREEDAPBCT) THEN
701: 701: 
702:                IF(NATOMS.NE.(NONEDAPBC**3)) THEN702:                IF(NATOMS.NE.(NONEDAPBC**3)) THEN
703:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'703:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'
704:                ENDIF704:                ENDIF
705: 705: 
706:                CALL ENERGY_3D_APBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)706:                CALL ENERGY_3D_APBC(COORDS,VNEW,ENERGY,GTEST,STEST)
707: 707: 
708: 708: 
709:                IF (PTEST) THEN709:                IF (PTEST) THEN
710:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '710:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
711:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '711:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
712:                ENDIF712:                ENDIF
713: 713: 
714:             ELSE IF( FOURDAPBCT) THEN714:             ELSE IF( FOURDAPBCT) THEN
715: 715: 
716:                IF(NATOMS.NE.(NONEDAPBC**4)) THEN716:                IF(NATOMS.NE.(NONEDAPBC**4)) THEN
717:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'717:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'
718:                ENDIF718:                ENDIF
719: 719: 
720:                CALL ENERGY_4D_APBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)720:                CALL ENERGY_4D_APBC(COORDS,VNEW,ENERGY,GTEST,STEST)
721: 721: 
722: 722: 
723:                IF (PTEST) THEN723:                IF (PTEST) THEN
724:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '724:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
725:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '725:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
726:                ENDIF726:                ENDIF
727: 727: 
728:             ELSE IF(FOURDPBCT) THEN728:             ELSE IF(FOURDPBCT) THEN
729: 729: 
730:                IF(NATOMS.NE.(NONEDAPBC**4)) THEN730:                IF(NATOMS.NE.(NONEDAPBC**4)) THEN
731:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'731:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number uunder variables in odata'
732:                ENDIF732:                ENDIF
733: 733: 
734:                CALL ENERGY_4D_PBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)734:                CALL ENERGY_4D_PBC(COORDS,VNEW,ENERGY,GTEST,STEST)
735: 735: 
736: 736: 
737:                IF (PTEST) THEN737:                IF (PTEST) THEN
738:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '738:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
739:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '739:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
740:                ENDIF740:                ENDIF
741: 741: 
742: 742: 
743:             ELSE IF(THREEDPBCT) THEN743:             ELSE IF(THREEDPBCT) THEN
744: 744: 
745:                IF(NATOMS.NE.(NONEDAPBC**3)) THEN745:                IF(NATOMS.NE.(NONEDAPBC**3)) THEN
746:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number under variables in odata'746:                   PRINT *, 'potential.f > Number of lattice sites specified does not equal number under variables in odata'
747:                ENDIF747:                ENDIF
748: 748: 
749:                CALL ENERGY_3D_PBC(COORDS,VNEW,ENERGY,GTEST,SSTEST)749:                CALL ENERGY_3D_PBC(COORDS,VNEW,ENERGY,GTEST,STEST)
750: 750: 
751: 751: 
752:                IF (PTEST) THEN752:                IF (PTEST) THEN
753:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '753:                   WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' '
754:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '754:                   WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' '
755:                ENDIF755:                ENDIF
756: 756: 
757:             END IF757:             END IF
758: 758: 
759:             ! IF (RESTART) THEN759:             ! IF (RESTART) THEN
795:             ! 430-441, 2009.795:             ! 430-441, 2009.
796:             ! 796:             ! 
797:          ELSE IF (ZSYM(NATOMS).EQ.'SB')then797:          ELSE IF (ZSYM(NATOMS).EQ.'SB')then
798: 798: 
799: 799: 
800:           IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN800:           IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
801:                XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)801:                XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)
802:                CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS,XRIGIDCOORDS)802:                CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS,XRIGIDCOORDS)
803:             ENDIF803:             ENDIF
804: 804: 
805:          IF (SSTEST) PRINT '(A)',' potential> ERROR - calling SBM with SSTEST true'805:          IF (STEST) PRINT '(A)',' potential> ERROR - calling SBM with STEST true'
806:          CALL SBM(COORDS,NATOMS,VNEW,ENERGY,GTEST,SSTEST)806:          CALL SBM(COORDS,NATOMS,VNEW,ENERGY,GTEST,STEST)
807:          IF (PTEST) THEN807:          IF (PTEST) THEN
808:             WRITE(*,10) ' Energy for last cycle=',ENERGY808:             WRITE(*,10) ' Energy for last cycle=',ENERGY
809:             WRITE(ESTRING,10) ' Energy for last cycle=',ENERGY809:             WRITE(ESTRING,10) ' Energy for last cycle=',ENERGY
810:          ENDIF810:          ENDIF
811: 811: 
812:             IF ( RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN812:             IF ( RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN
813:                IF (SSTEST) THEN813:                IF (SSTEST) THEN
814:                   CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS,XRIGIDHESS, RBAANORMALMODET)814:                   CALL TRANSFORMHESSIAN(HESS, VNEW, XRIGIDCOORDS,XRIGIDHESS, RBAANORMALMODET)
815:                   HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0815:                   HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0
816:                   HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0816:                   HESS(:,DEGFREEDOMS+1:3*NATOMS) = 0.0D0
995:                WRITE(*,20) ' Three-body contribution= ',P3995:                WRITE(*,20) ' Three-body contribution= ',P3
996:                WRITE(*,20) ' Z parameter=             ',ZSTAR996:                WRITE(*,20) ' Z parameter=             ',ZSTAR
997:             ENDIF997:             ENDIF
998:          ELSE IF (ZSYM(NATOMS).EQ.'Z1') THEN998:          ELSE IF (ZSYM(NATOMS).EQ.'Z1') THEN
999:             CALL Z1(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST,PARAM1,PARAM2,PARAM3)999:             CALL Z1(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST,PARAM1,PARAM2,PARAM3)
1000:             IF (PTEST) THEN1000:             IF (PTEST) THEN
1001:                WRITE(*,20) ' potential> Energy for last cycle=',ENERGY1001:                WRITE(*,20) ' potential> Energy for last cycle=',ENERGY
1002:                WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY1002:                WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY
1003:             ENDIF1003:             ENDIF
1004:          ELSE IF (ZSYM(NATOMS).EQ.'ZF') THEN1004:          ELSE IF (ZSYM(NATOMS).EQ.'ZF') THEN
1005:             CALL Z2FASTER(NATOMS,COORDS,PARAM1,PARAM2,PARAM3,ENERGY,VNEW,SSTEST)1005:             CALL Z2FASTER(NATOMS,COORDS,PARAM1,PARAM2,PARAM3,ENERGY,VNEW,STEST)
1006:             IF (PTEST) THEN1006:             IF (PTEST) THEN
1007:                WRITE(*,20) ' potential> Energy for last cycle=',ENERGY1007:                WRITE(*,20) ' potential> Energy for last cycle=',ENERGY
1008:                WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY1008:                WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY
1009:             ENDIF1009:             ENDIF
1010:          ELSE IF (ZSYM(NATOMS).EQ.'Z2') THEN1010:          ELSE IF (ZSYM(NATOMS).EQ.'Z2') THEN
1011:             CALL Z2(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST,PARAM1,PARAM2,PARAM3)1011:             CALL Z2(NATOMS, COORDS, VNEW, ENERGY, GTEST, SSTEST,PARAM1,PARAM2,PARAM3)
1012:             IF (PTEST) THEN1012:             IF (PTEST) THEN
1013:                WRITE(*,20) ' potential> Energy for last cycle=',ENERGY1013:                WRITE(*,20) ' potential> Energy for last cycle=',ENERGY
1014:                WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY1014:                WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY
1015:             ENDIF1015:             ENDIF
3203: ! AMBER 9 energy calculation3203: ! AMBER 9 energy calculation
3204:          ELSE IF (AMBERT) THEN3204:          ELSE IF (AMBERT) THEN
3205:             VNEW = 0.0D03205:             VNEW = 0.0D0
3206:             IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN3206:             IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
3207:                XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)3207:                XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)
3208:                CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)3208:                CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)
3209:             ENDIF3209:             ENDIF
3210:             IF (CHECKCISTRANSALWAYS) CALL CHECK_CISTRANS_PROTEIN(COORDS,NATOMS,goodstructure1,MINOMEGA,CISARRAY1)3210:             IF (CHECKCISTRANSALWAYS) CALL CHECK_CISTRANS_PROTEIN(COORDS,NATOMS,goodstructure1,MINOMEGA,CISARRAY1)
3211:             IF (CHECKCISTRANSALWAYSDNA) CALL CHECK_CISTRANS_DNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)3211:             IF (CHECKCISTRANSALWAYSDNA) CALL CHECK_CISTRANS_DNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)
3212:             IF (CHECKCISTRANSALWAYSRNA) CALL CHECK_CISTRANS_RNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)3212:             IF (CHECKCISTRANSALWAYSRNA) CALL CHECK_CISTRANS_RNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)
3213:             CALL AMBERENERGIES(COORDS,GRADATOMS,ENERGY,GTEST,SSTEST)3213:             CALL AMBERENERGIES(COORDS,GRADATOMS,ENERGY,GTEST,STEST)
3214:             VNEW(1:3*NATOMS) = GRADATOMS(:)3214:             VNEW(1:3*NATOMS) = GRADATOMS(:)
3215:             IF (STEST) CALL AMBERSECDER(COORDS,.TRUE.)3215:             IF (STEST) CALL AMBERSECDER(COORDS,.TRUE.)
3216:             IF (PTEST) THEN3216:             IF (PTEST) THEN
3217:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' kcal/mol'3217:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' kcal/mol'
3218:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' kcal/mol'3218:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' kcal/mol'
3219:             ENDIF3219:             ENDIF
3220:             IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT) ) THEN3220:             IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT) ) THEN
3221:                IF (STEST) THEN3221:                IF (STEST) THEN
3222:                   CALL TRANSFORMHESSIAN(HESS, GRADATOMS, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)3222:                   CALL TRANSFORMHESSIAN(HESS, GRADATOMS, XRIGIDCOORDS, XRIGIDHESS, RBAANORMALMODET)
3223:                   HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D03223:                   HESS(DEGFREEDOMS+1:3*NATOMS,:) = 0.0D0
3240:             ! NAB structures initialised in keywords.f3240:             ! NAB structures initialised in keywords.f
3241:             VNEW = 0.0D03241:             VNEW = 0.0D0
3242:             IF(STEST) THEN3242:             IF(STEST) THEN
3243:                ! Analytical second derivatives3243:                ! Analytical second derivatives
3244:                ! hk2863244:                ! hk286
3245:                IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN3245:                IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
3246:                   XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)3246:                   XRIGIDCOORDS(1:DEGFREEDOMS) = COORDS(1:DEGFREEDOMS)
3247:                   CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)3247:                   CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)
3248:                ENDIF3248:                ENDIF
3249:                ! hk286 - numerical second derivatives3249:                ! hk286 - numerical second derivatives
3250:                ! CALL NABSECDER(COORDS,SSTEST)3250:                ! CALL NABSECDER(COORDS,STEST)
3251: 3251: 
3252:                IF (ALLOCATED(HESS)) DEALLOCATE(HESS)3252:                IF (ALLOCATED(HESS)) DEALLOCATE(HESS)
3253:                IF (.NOT.ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))3253:                IF (.NOT.ALLOCATED(TEMPHESS)) ALLOCATE(TEMPHESS(9*NATOMS*NATOMS))
3254:                TEMPHESS(:) = 0.0D03254:                TEMPHESS(:) = 0.0D0
3255:                CALL MME2WRAPPER(COORDS,ENERGY,GRADATOMS,TEMPHESS,ATMASS,GRAD1)3255:                CALL MME2WRAPPER(COORDS,ENERGY,GRADATOMS,TEMPHESS,ATMASS,GRAD1)
3256:                IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN3256:                IF (SQRT(ABS(ENERGY)) > 1.0D40) THEN
3257:                   WRITE(*, *) "Linear dihedral detected in NAB routines (sff2.c)."3257:                   WRITE(*, *) "Linear dihedral detected in NAB routines (sff2.c)."
3258:                END IF3258:                END IF
3259:                VNEW(1:3*NATOMS) = GRADATOMS(:)3259:                VNEW(1:3*NATOMS) = GRADATOMS(:)
3260:                ALLOCATE(HESS(3*NATOMS,3*NATOMS))3260:                ALLOCATE(HESS(3*NATOMS,3*NATOMS))
3285:                   CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)3285:                   CALL TRANSFORMRIGIDTOC(1, NRIGIDBODY, COORDS, XRIGIDCOORDS)
3286:                ENDIF3286:                ENDIF
3287:                ! CALL AMBERNUMFIRSTDER(COORDS,GTEST)3287:                ! CALL AMBERNUMFIRSTDER(COORDS,GTEST)
3288: 3288: 
3289:                ! check cis-trans isomerisation for DNA or RNA3289:                ! check cis-trans isomerisation for DNA or RNA
3290:                IF (CHECKCISTRANSALWAYSDNA) CALL CHECK_CISTRANS_DNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)3290:                IF (CHECKCISTRANSALWAYSDNA) CALL CHECK_CISTRANS_DNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)
3291:                IF (CHECKCISTRANSALWAYSRNA) CALL CHECK_CISTRANS_RNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)3291:                IF (CHECKCISTRANSALWAYSRNA) CALL CHECK_CISTRANS_RNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)
3292:                IF (CHECKCISTRANSALWAYS) CALL CHECK_CISTRANS_PROTEIN(COORDS,NATOMS,goodstructure1,MINOMEGA,CISARRAY1)3292:                IF (CHECKCISTRANSALWAYS) CALL CHECK_CISTRANS_PROTEIN(COORDS,NATOMS,goodstructure1,MINOMEGA,CISARRAY1)
3293:                ! CALL CHECK_CISTRANS_DNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)3293:                ! CALL CHECK_CISTRANS_DNA(COORDS,NATOMS,ZSYM,GOODSTRUCTURE1)
3294:                ! STOP3294:                ! STOP
3295:                CALL AMBERENERGIES(COORDS,GRADATOMS,ENERGY,GTEST,SSTEST)3295:                CALL AMBERENERGIES(COORDS,GRADATOMS,ENERGY,GTEST,STEST)
3296:                VNEW(1:3*NATOMS) = GRADATOMS(:)3296:                VNEW(1:3*NATOMS) = GRADATOMS(:)
3297:                ! CALL MME(ENERGY,COORDS,VNEW,1)3297:                ! CALL MME(ENERGY,COORDS,VNEW,1)
3298:             END IF3298:             END IF
3299:             IF (PTEST) THEN3299:             IF (PTEST) THEN
3300:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' kcal/mol'3300:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' kcal/mol'
3301:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' kcal/mol'3301:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' kcal/mol'
3302:             ENDIF3302:             ENDIF
3303:             ! hk2863303:             ! hk286
3304:             IF (.NOT. ATOMRIGIDCOORDT) THEN3304:             IF (.NOT. ATOMRIGIDCOORDT) THEN
3305:                CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)3305:                CALL TRANSFORMGRAD(GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD)
3747:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' epsilon'3747:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,' epsilon'
3748:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' epsilon'3748:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,' epsilon'
3749:             ENDIF3749:             ENDIF
3750:          ENDIF3750:          ENDIF
3751:          IF (AAORIENTT) THEN3751:          IF (AAORIENTT) THEN
3752:             ! IF (DEBUG) PRINT '(A)',' potential> Adding additional angle-axis potential'3752:             ! IF (DEBUG) PRINT '(A)',' potential> Adding additional angle-axis potential'
3753:             ! ENERGY=0.0D0  !!!!!!!!!!!!!!!!!!!!! debug3753:             ! ENERGY=0.0D0  !!!!!!!!!!!!!!!!!!!!! debug
3754:             ! VNEW(1:NOPT)=0.0D0  !!!!!!!!!!!!!!!!!!!!! debug3754:             ! VNEW(1:NOPT)=0.0D0  !!!!!!!!!!!!!!!!!!!!! debug
3755:             ! HESS(1:NOPT,1:NOPT)=0.0D0  !!!!!!!!!!!!!!!!!!!!! debug3755:             ! HESS(1:NOPT,1:NOPT)=0.0D0  !!!!!!!!!!!!!!!!!!!!! debug
3756:             IF (SIGMAAA .EQ. 0.0D0) THEN3756:             IF (SIGMAAA .EQ. 0.0D0) THEN
3757:                CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,SSTEST)3757:                CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,STEST)
3758:             ELSE3758:             ELSE
3759:                CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,SSTEST)3759:                CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,STEST)
3760:             END IF3760:             END IF
3761:             IF (PTEST) THEN3761:             IF (PTEST) THEN
3762:                WRITE(*,10) ' potential> Energy for last cycle with additional angle-axis potential=',ENERGY3762:                WRITE(*,10) ' potential> Energy for last cycle with additional angle-axis potential=',ENERGY
3763:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY3763:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY
3764:             ENDIF3764:             ENDIF
3765:             ! Debug tools3765:             ! Debug tools
3766:             ! 3766:             ! 
3767:             ! DIFF=1.0D-43767:             ! DIFF=1.0D-4
3768:             ! PRINT*,'analytic and numerical gradients: NATOMS=',NATOMS3768:             ! PRINT*,'analytic and numerical gradients: NATOMS=',NATOMS
3769:             ! DO J1=1,3*NATOMS3769:             ! DO J1=1,3*NATOMS
3770:             ! EPLUS=0.0D03770:             ! EPLUS=0.0D0
3771:             ! EMINUS=0.0D03771:             ! EMINUS=0.0D0
3772:             ! VPLUS(1:NOPT)=0.0D03772:             ! VPLUS(1:NOPT)=0.0D0
3773:             ! VMINUS(1:NOPT)=0.0D03773:             ! VMINUS(1:NOPT)=0.0D0
3774:             ! COORDS(J1)=COORDS(J1)+DIFF3774:             ! COORDS(J1)=COORDS(J1)+DIFF
3775:             ! IF (SIGMAAA .EQ. 0.0D0) THEN3775:             ! IF (SIGMAAA .EQ. 0.0D0) THEN
3776:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,SSTEST)3776:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,STEST)
3777:             ! ELSE3777:             ! ELSE
3778:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,SSTEST)3778:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,STEST)
3779:             ! END IF3779:             ! END IF
3780:             ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF3780:             ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF
3781:             ! IF (SIGMAAA .EQ. 0.0D0) THEN3781:             ! IF (SIGMAAA .EQ. 0.0D0) THEN
3782:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,SSTEST)3782:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,STEST)
3783:             ! ELSE3783:             ! ELSE
3784:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,SSTEST)3784:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,STEST)
3785:             ! END IF3785:             ! END IF
3786:             ! COORDS(J1)=COORDS(J1)+DIFF3786:             ! COORDS(J1)=COORDS(J1)+DIFF
3787:             ! IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.0.0D0)) THEN3787:             ! IF ((ABS(VNEW(J1)).NE.0.0D0).AND.(100.0D0*ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.0.0D0)) THEN
3788:             ! IF (ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.1.0D-2) THEN3788:             ! IF (ABS((VNEW(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/VNEW(J1)).GT.1.0D-2) THEN
3789:             ! WRITE(*,'(I5,3G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF),(EPLUS-EMINUS)/(2.0D0*DIFF*VNEW(J1)),'   X'3789:             ! WRITE(*,'(I5,3G20.10,A)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF),(EPLUS-EMINUS)/(2.0D0*DIFF*VNEW(J1)),'   X'
3790:             ! ELSE3790:             ! ELSE
3791:             ! WRITE(*,'(I5,3G20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF),(EPLUS-EMINUS)/(2.0D0*DIFF*VNEW(J1))3791:             ! WRITE(*,'(I5,3G20.10)') J1,VNEW(J1),(EPLUS-EMINUS)/(2.0D0*DIFF),(EPLUS-EMINUS)/(2.0D0*DIFF*VNEW(J1))
3792:             ! ENDIF3792:             ! ENDIF
3793:             ! ENDIF3793:             ! ENDIF
3794:             ! ENDDO3794:             ! ENDDO
3795:             ! PRINT*,'analytic and numerical second derivatives:'3795:             ! PRINT*,'analytic and numerical second derivatives:'
3796:             ! DO J1=1,3*NATOMS3796:             ! DO J1=1,3*NATOMS
3797:             ! EPLUS=0.0D03797:             ! EPLUS=0.0D0
3798:             ! EMINUS=0.0D03798:             ! EMINUS=0.0D0
3799:             ! VPLUS(1:NOPT)=0.0D03799:             ! VPLUS(1:NOPT)=0.0D0
3800:             ! VMINUS(1:NOPT)=0.0D03800:             ! VMINUS(1:NOPT)=0.0D0
3801:             ! COORDS(J1)=COORDS(J1)+DIFF3801:             ! COORDS(J1)=COORDS(J1)+DIFF
3802:             ! IF (SIGMAAA .EQ. 0.0D0) THEN3802:             ! IF (SIGMAAA .EQ. 0.0D0) THEN
3803:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,SSTEST)3803:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,STEST)
3804:             ! ELSE3804:             ! ELSE
3805:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,SSTEST)3805:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,STEST)
3806:             ! END IF3806:             ! END IF
3807:             ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF3807:             ! COORDS(J1)=COORDS(J1)-2.0D0*DIFF
3808:             ! IF (SIGMAAA .EQ. 0.0D0) THEN3808:             ! IF (SIGMAAA .EQ. 0.0D0) THEN
3809:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,SSTEST)3809:             ! CALL AAORIENT(COORDS,VNEW,ENERGY,GTEST,STEST)
3810:             ! ELSE3810:             ! ELSE
3811:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,SSTEST)3811:             ! CALL AAORIENTSR(COORDS,VNEW,ENERGY,GTEST,STEST)
3812:             ! END IF3812:             ! END IF
3813:             ! COORDS(J1)=COORDS(J1)+DIFF3813:             ! COORDS(J1)=COORDS(J1)+DIFF
3814:             ! DO J2=1,3*NATOMS3814:             ! DO J2=1,3*NATOMS
3815:             ! IF (ABS(HESS(J1,J2)).GT.0.0D0) THEN3815:             ! IF (ABS(HESS(J1,J2)).GT.0.0D0) THEN
3816:             ! IF (ABS((HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D-2) THEN3816:             ! IF (ABS((HESS(J1,J2)-(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF))/HESS(J1,J2)).GT.1.0D-2) THEN
3817:             ! WRITE(*,'(2I5,3F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),3817:             ! WRITE(*,'(2I5,3F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),
3818:             ! &                                                        ABS((VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF*HESS(J1,J2))),'   X'3818:             ! &                                                        ABS((VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF*HESS(J1,J2))),'   X'
3819:             ! ELSE3819:             ! ELSE
3820:             ! WRITE(*,'(2I5,3F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),3820:             ! WRITE(*,'(2I5,3F20.10,A)') J1,J2,HESS(J1,J2),(VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF),
3821:             ! &                                                        ABS((VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF*HESS(J1,J2)))3821:             ! &                                                        ABS((VPLUS(J2)-VMINUS(J2))/(2.0D0*DIFF*HESS(J1,J2)))
3924:          ! WRITE(*,20) ' potential> Energy for last cycle=',ENERGY3924:          ! WRITE(*,20) ' potential> Energy for last cycle=',ENERGY
3925:          ! WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY3925:          ! WRITE(ESTRING,20) ' potential> Energy for last cycle=',ENERGY
3926:          ! WRITE(*,'(A23,7X,2F20.10)') ' RHO and DELTA=',PARAM1, PARAM23926:          ! WRITE(*,'(A23,7X,2F20.10)') ' RHO and DELTA=',PARAM1, PARAM2
3927:          ! ENDIF3927:          ! ENDIF
3928:          ! ENDIF3928:          ! ENDIF
3929:          ! 3929:          ! 
3930:          ! Double well potential between the first two atoms.3930:          ! Double well potential between the first two atoms.
3931:          ! 3931:          ! 
3932:          IF (DOUBLET) THEN3932:          IF (DOUBLET) THEN
3933:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'3933:             PRINT*,'WARNING - this potential has not been tested in OPTIM.3.0'
3934:             CALL DOUBLE(NATOMS,COORDS,VNEW,EDOUBLE,GTEST,SSTEST,PARAM4,PARAM5,PARAM6)3934:             CALL DOUBLE(NATOMS,COORDS,VNEW,EDOUBLE,GTEST,STEST,PARAM4,PARAM5,PARAM6)
3935:             ENERGY=ENERGY+EDOUBLE3935:             ENERGY=ENERGY+EDOUBLE
3936:             IF (PTEST) THEN3936:             IF (PTEST) THEN
3937:                WRITE(*,'(A,F20.10,A)') ' potential> Energy for last cycle including double well=     ',ENERGY,' epsilon'3937:                WRITE(*,'(A,F20.10,A)') ' potential> Energy for last cycle including double well=     ',ENERGY,' epsilon'
3938:                WRITE(ESTRING,'(A,F20.10,A)') ' potential> Energy for last cycle including double well=     ',ENERGY,' epsilon'3938:                WRITE(ESTRING,'(A,F20.10,A)') ' potential> Energy for last cycle including double well=     ',ENERGY,' epsilon'
3939:             ENDIF3939:             ENDIF
3940:          ENDIF3940:          ENDIF
3941: 3941: 
3942:          IF (GAUSSIAN.OR.CADPAC.OR.GAMESSUK.OR.GAMESSUS) THEN3942:          IF (GAUSSIAN.OR.CADPAC.OR.GAMESSUK.OR.GAMESSUS) THEN
3943:             INQUIRE(FILE='abenergy',EXIST=ETEST)3943:             INQUIRE(FILE='abenergy',EXIST=ETEST)
3944:             IF (ETEST) THEN3944:             IF (ETEST) THEN
4090:             IF (GTEST) VNEW(1:NOPT)=-VNEW(1:NOPT)4090:             IF (GTEST) VNEW(1:NOPT)=-VNEW(1:NOPT)
4091:             IF (SSTEST) HESS(1:NOPT,1:NOPT)=-HESS(1:NOPT,1:NOPT)4091:             IF (SSTEST) HESS(1:NOPT,1:NOPT)=-HESS(1:NOPT,1:NOPT)
4092:          ENDIF4092:          ENDIF
4093: 4093: 
4094:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN4094:       IF (RIGIDINIT .AND. (ATOMRIGIDCOORDT .EQV. .FALSE.) ) THEN
4095:         ! sn402: transform back to rigid coords if required4095:         ! sn402: transform back to rigid coords if required
4096:           COORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)4096:           COORDS(1:DEGFREEDOMS) = XRIGIDCOORDS(1:DEGFREEDOMS)
4097:           COORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D04097:           COORDS(DEGFREEDOMS+1:3*NATOMS) = 0.0D0
4098:       ENDIF4098:       ENDIF
4099: 4099: 
4100: !          WRITE(*,'(A,G30.20)') 'Energy in potential:',ENERGY4100: !          WRITE(*,'(A,G30.20)') 'Energy in POTENTIAL:',ENERGY
4101: !          PRINT*,'GTEST,STEST,SSTEST=',GTEST,STEST,SSTEST4101: !          PRINT*,'GTEST,SSTEST=',GTEST,SSTEST
4102: !          WRITE(*,'(A,G20.10)') 'RMS in potential=',RMS4102: !          WRITE(*,'(A,G20.10)') 'RMS in potential=',RMS
4103: !          PRINT*,'coords in potential:'4103: !          PRINT*,'coords in potential:'
4104: !          WRITE(*,'(6G25.15)') (COORDS(J1),J1=1,NOPT)4104: !          WRITE(*,'(6G25.15)') (COORDS(J1),J1=1,NOPT)
4105:          ! CALL FLUSH(6)4105:          ! CALL FLUSH(6)
4106:          ! PRINT*,'PARAMS'4106:          ! PRINT*,'PARAMS'
4107:          ! WRITE(*,'(3F20.10)') PARAM1,PARAM2,PARAM34107:          ! WRITE(*,'(3F20.10)') PARAM1,PARAM2,PARAM3
4108: !          IF (GTEST) PRINT*,'grad:'4108: !          IF (GTEST) PRINT*,'grad:'
4109: !          IF (GTEST) WRITE(*,'(6F15.5)') (VNEW(J1),J1=1,NOPT)4109: !          IF (GTEST) WRITE(*,'(6F15.5)') (VNEW(J1),J1=1,NOPT)
4110: !          IF (SSTEST) PRINT*,'hess:'4110:          ! IF (SSTEST) PRINT*,'hess:'
4111: !          IF (SSTEST) WRITE(*,'(6F15.5)') ((HESS(J1,J2),J1=1,NOPT),J2=1,NOPT)4111:          ! IF (SSTEST) WRITE(*,'(6F15.5)') ((HESS(J1,J2),J1=1,NOPT),J2=1,NOPT)
4112: 4112: 
4113:          CALL MYCPU_TIME(TIME,.FALSE.)4113:          CALL MYCPU_TIME(TIME,.FALSE.)
4114:          IF (SSTEST) THEN4114:          IF (SSTEST) THEN
4115:             SCALL=SCALL+14115:             SCALL=SCALL+1
4116:             STIME=STIME+TIME-TIME04116:             STIME=STIME+TIME-TIME0
4117:          ELSE IF (GTEST) THEN4117:          ELSE IF (GTEST) THEN
4118:             FCALL=FCALL+14118:             FCALL=FCALL+1
4119:             FTIME=FTIME+TIME-TIME04119:             FTIME=FTIME+TIME-TIME0
4120:          ELSE4120:          ELSE
4121:             ECALL=ECALL+14121:             ECALL=ECALL+1


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0