hdiff output

r33407/congrad.f90 2017-10-23 10:30:21.029307004 +0100 r33406/congrad.f90 2017-10-23 10:30:22.717329338 +0100
508: !508: !
509: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)509: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
510: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &510: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
511:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &511:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
512:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &512:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &
513:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &513:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &
514:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET514:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET
515: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG515: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
516: IMPLICIT NONE516: IMPLICIT NONE
517:            517:            
518: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2),JMAX,IMAX518: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2)
519: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX519: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS
520: INTEGER JJMAX(INTIMAGE+2) 
521: DOUBLE PRECISION  EEMAX(INTIMAGE+2) 
522: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1520: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
523: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)521: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
524: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI522: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
525: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)523: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)
526: LOGICAL NOINT, LPRINT524: LOGICAL NOINT, LPRINT
527: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)525: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)
528: LOGICAL IMGFREEZE(INTIMAGE), PRINTE526: LOGICAL IMGFREEZE(INTIMAGE), PRINTE
529: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), CCLOCAL527: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), CCLOCAL
530: 528: 
531: PRINTE=.FALSE.529: PRINTE=.FALSE.
545: !543: !
546: !  Constraint energy and forces.544: !  Constraint energy and forces.
547: !545: !
548: ! For J1 we consider the line segment between image J1-1 and J1.546: ! For J1 we consider the line segment between image J1-1 and J1.
549: ! There are INTIMAGE+1 line segments in total, with an energy contribution547: ! There are INTIMAGE+1 line segments in total, with an energy contribution
550: ! and corresponding gradient terms for each. 548: ! and corresponding gradient terms for each. 
551: ! A and B refer to atoms, 1 and 2 to images J1-1 and J1 corresponding to J1-2 and J1-1 below.549: ! A and B refer to atoms, 1 and 2 to images J1-1 and J1 corresponding to J1-2 and J1-1 below.
552: !550: !
553: ! IMGFREEZE(1:INTIMAGE) refers to the images excluding end points!551: ! IMGFREEZE(1:INTIMAGE) refers to the images excluding end points!
554: !552: !
555: EMAX=-1.0D200 
556: EEMAX(1:INTIMAGE+2)=-1.0D200 
557: JJMAX(1:INTIMAGE+2)=-1 
558: JMAX=-1 
559: IMAX=-1 
560: DO J2=1,NCONSTRAINT553: DO J2=1,NCONSTRAINT
561:    IF (.NOT.CONACTIVE(J2)) CYCLE554:    IF (.NOT.CONACTIVE(J2)) CYCLE
562:    CCLOCAL=CONCUTLOCAL(J2)555:    CCLOCAL=CONCUTLOCAL(J2)
563:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS556:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
564:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)557:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
565: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG558: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG
566:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN559:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN
567:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)560:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)
568:       STOP561:       STOP
569:    ENDIF562:    ENDIF
604: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)597: !     CONCUT=CONCUTFRAC*CONDISTREF(J2)
605: !598: !
606: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.599: ! terms for image J1 - non-zero derivatives only for J1. D2 is the distance for image J1.
607: !600: !
608:       IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &601:       IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &
609:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL602:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL
610:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 603:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 
611:          DUMMY=D2-CONDISTREFLOCAL(J2)  604:          DUMMY=D2-CONDISTREFLOCAL(J2)  
612:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)605:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
613:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)606:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
614:          IF (DUMMY.GT.EMAX) THEN 
615:             IMAX=J1 
616:             JMAX=J2 
617:             EMAX=DUMMY 
618:          ENDIF 
619:          IF (DUMMY.GT.EEMAX(J1)) THEN 
620:             JJMAX(J1)=J2 
621:             EEMAX(J1)=DUMMY 
622:          ENDIF 
623:          EEE(J1)=EEE(J1)+DUMMY607:          EEE(J1)=EEE(J1)+DUMMY
624:          CONE(J1)=CONE(J1)+DUMMY608:          CONE(J1)=CONE(J1)+DUMMY
625:          ECON=ECON      +DUMMY609:          ECON=ECON      +DUMMY
626:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &610:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
627:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)611:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
628:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)612:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
629:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)613:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
630:       ENDIF614:       ENDIF
631: !615: !
632: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.616: ! Don't add energy contributions to EEE(2) from D1, since the gradients are non-zero only for image 1.
633: !617: !
634: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.618: ! terms for image J1-1 - non-zero derivatives only for J1-1. D1 is the distance for image J1-1.
635: !619: !
636: !     IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &620: !     IF (LPRINT) WRITE(*, '(A,I6,5G15.5)') &
637: ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL621: ! &       'J1,D2,D1,DINT,MAX diff,CCLOCAL=',J1,D2,D1,DINT,ABS(D1-CONDISTREFLOCAL(J2)),CCLOCAL
638:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  622:       IF ((ABS(D1-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.GT.2)) THEN  
639:          DUMMY=D1-CONDISTREFLOCAL(J2)  623:          DUMMY=D1-CONDISTREFLOCAL(J2)  
640:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)624:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1(1:3)
641:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)625:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
642: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN626:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
643: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &627:             WRITE(*, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
644: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY628:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
645: !        ENDIF 
646:          IF (DUMMY.GT.EMAX) THEN 
647:             IMAX=J1 
648:             JMAX=J2 
649:             EMAX=DUMMY 
650:          ENDIF 
651:          IF (DUMMY.GT.EEMAX(J1-1)) THEN 
652:             JJMAX(J1-1)=J2 
653:             EEMAX(J1-1)=DUMMY 
654:          ENDIF629:          ENDIF
655:          EEE(J1-1)=EEE(J1-1)+DUMMY630:          EEE(J1-1)=EEE(J1-1)+DUMMY
656:          CONE(J1-1)=CONE(J1-1)+DUMMY631:          CONE(J1-1)=CONE(J1-1)+DUMMY
657:          ECON=ECON      +DUMMY632:          ECON=ECON      +DUMMY
658:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &633:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
659:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)634:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
660:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)635:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
661:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)636:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
662:       ENDIF637:       ENDIF
663:       IF ((.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN638:       IF ((.NOT.NOINT).AND.(ABS(DINT-CONDISTREFLOCAL(J2)).GT.CCLOCAL)) THEN
664:          DUMMY=DINT-CONDISTREFLOCAL(J2)  639:          DUMMY=DINT-CONDISTREFLOCAL(J2)  
665:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)640:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G1INT(1:3)
666:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)641:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
667:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)642:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
668:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)643:          REPGRAD(1:3)=2*INTMINFAC*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2INT(1:3)
669:          DUMMY=INTMINFAC*INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)644:          DUMMY=INTMINFAC*INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
670: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN645:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
671: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &646:             WRITE(*, '(A,2I6,2L5,G20.10)') 'B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
672: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY647:   &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
673: !        ENDIF 
674:          ECON=ECON+DUMMY 
675:          IF (DUMMY.GT.EMAX) THEN 
676:             IMAX=J1 
677:             JMAX=J2 
678:             EMAX=DUMMY 
679:          ENDIF 
680:          IF (DUMMY.GT.EEMAX(J1-1)) THEN 
681:             JJMAX(J1-1)=J2 
682:             EEMAX(J1-1)=DUMMY 
683:          ENDIF 
684:          IF (DUMMY.GT.EEMAX(J1)) THEN 
685:             JJMAX(J1)=J2 
686:             EEMAX(J1)=DUMMY 
687:          ENDIF648:          ENDIF
 649:          ECON=ECON+DUMMY
688:          IF (J1.EQ.2) THEN650:          IF (J1.EQ.2) THEN
689:             EEE(J1)=EEE(J1)+DUMMY651:             EEE(J1)=EEE(J1)+DUMMY
690:             CONEINT(J1)=CONEINT(J1)+DUMMY652:             CONEINT(J1)=CONEINT(J1)+DUMMY
691:             NCONINT(J1)=NCONINT(J1)+1653:             NCONINT(J1)=NCONINT(J1)+1
692:          ELSE IF (J1.LT.INTIMAGE+2) THEN654:          ELSE IF (J1.LT.INTIMAGE+2) THEN
693:             EEE(J1)=EEE(J1)+DUMMY/2.0D0655:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
694:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0656:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
695:             CONEINT(J1)=CONEINT(J1)+DUMMY/2.0D0657:             CONEINT(J1)=CONEINT(J1)+DUMMY/2.0D0
696:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY/2.0D0658:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY/2.0D0
697:             NCONINT(J1)=NCONINT(J1)+1659:             NCONINT(J1)=NCONINT(J1)+1
701:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY663:             CONEINT(J1-1)=CONEINT(J1-1)+DUMMY
702:             NCONINT(J1-1)=NCONINT(J1-1)+1664:             NCONINT(J1-1)=NCONINT(J1-1)+1
703:          ENDIF665:          ENDIF
704: !        WRITE(*, '(A,4I6,G15.5)') 'in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &666: !        WRITE(*, '(A,4I6,G15.5)') 'in2 J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
705: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)667: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
706:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)668:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
707:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)669:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
708:       ENDIF670:       ENDIF
709:    ENDDO671:    ENDDO
710: ENDDO672: ENDDO
711: ! DO J2=1,INTIMAGE+2 
712: !    IF (JJMAX(J2).GT.0) THEN 
713: !       WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest constraint contribution for image ',J2,' constraint ',JJMAX(J2),' atoms ', & 
714: !  &                                CONI(JJMAX(J2)),CONJ(JJMAX(J2)) 
715: !    ENDIF 
716: ! ENDDO 
717: IF (JMAX.GT.0) THEN 
718:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest constraint contribution for any image in image ',IMAX, & 
719:  & ' constraint ',JMAX, & 
720:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX) 
721: ENDIF 
722: 673: 
723: ! INTCONST=INTCONSTRAINREPCUT**13674: ! INTCONST=INTCONSTRAINREPCUT**13
724: 675: 
725: EMAX=-1.0D200 
726: EEMAX(1:INTIMAGE+2)=-1.0D200 
727: JJMAX(1:INTIMAGE+2)=-1 
728: JMAX=-1 
729: IMAX=-1 
730: DO J2=1,NNREPULSIVE676: DO J2=1,NNREPULSIVE
731: !  INTCONST=NREPCUT(J2)**13677: !  INTCONST=NREPCUT(J2)**13
732:    INTCONST=NREPCUT(J2)**3678:    INTCONST=NREPCUT(J2)**3
733:    DO J1=2,INTIMAGE+2679:    DO J1=2,INTIMAGE+2
734:       IF (FREEZENODEST) THEN680:       IF (FREEZENODEST) THEN
735:          IF (J1.EQ.2) THEN681:          IF (J1.EQ.2) THEN
736:             IF (IMGFREEZE(1)) CYCLE682:             IF (IMGFREEZE(1)) CYCLE
737:          ELSE IF (J1.EQ.INTIMAGE+2) THEN683:          ELSE IF (J1.EQ.INTIMAGE+2) THEN
738:             IF (IMGFREEZE(INTIMAGE)) CYCLE684:             IF (IMGFREEZE(INTIMAGE)) CYCLE
739:          ELSE685:          ELSE
740:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE686:             IF (IMGFREEZE(J1-2).AND.IMGFREEZE(J1-1)) CYCLE
741:          ENDIF687:          ENDIF
742:       ENDIF688:       ENDIF
743:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN689:       IF (INTFROZEN(NREPI(J2)).AND.INTFROZEN(NREPJ(J2))) THEN
744:          WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)690:          WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** repulsion ',J2,' between frozen atoms ',NREPI(J2),NREPJ(J2)
745:          STOP691:          STOP
746:       ENDIF692:       ENDIF
 693: !     WRITE(*,'(A,2I8,6G20.10)') 'congrad2> B J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)
747:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)694:       NI1=(3*NATOMS)*(J1-2)+3*(NREPI(J2)-1)
748:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)695:       NI2=(3*NATOMS)*(J1-1)+3*(NREPI(J2)-1)
749:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)696:       NJ1=(3*NATOMS)*(J1-2)+3*(NREPJ(J2)-1)
750:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)697:       NJ2=(3*NATOMS)*(J1-1)+3*(NREPJ(J2)-1)
751:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)698:       R1AX=XYZ(NI1+1); R1AY=XYZ(NI1+2); R1AZ=XYZ(NI1+3)
752:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)699:       R1BX=XYZ(NJ1+1); R1BY=XYZ(NJ1+2); R1BZ=XYZ(NJ1+3)
753:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)700:       R2AX=XYZ(NI2+1); R2AY=XYZ(NI2+2); R2AZ=XYZ(NI2+3)
754:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)701:       R2BX=XYZ(NJ2+1); R2BY=XYZ(NJ2+2); R2BZ=XYZ(NJ2+3)
755: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN702: !     IF (r2ax**2+r2ay**2+r2az**2+r2bx**2+r2by**2+r2bz**2-2*(r2ax*r2bx+r2ay*r2by+r2az*r2bz).EQ.0.0D0) THEN
756:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-50) THEN703:       IF ((r1ax-r1bx-r2ax+r2bx)**2+(r1ay-r1by-r2ay+r2by)**2+(r1az-r1bz-r2az+r2bz)**2.LT.1.0D-50) THEN
761: !        WRITE(*,'(A,7I10)') 'frames ',J1-1,J1708: !        WRITE(*,'(A,7I10)') 'frames ',J1-1,J1
762:          D1=1.0D100; D2=1.0D100; NOINT=.TRUE.  ! to skip the next blocks709:          D1=1.0D100; D2=1.0D100; NOINT=.TRUE.  ! to skip the next blocks
763:       ELSE710:       ELSE
764:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &711:          CALL MINMAXD2R(R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ, &
765:   &                    D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))712:   &                    D2,D1,DINT,DSQ2,DSQ1,DSQI,G1,G2,G1INT,G2INT,NOINT,.FALSE.,NREPCUT(J2))
766:       ENDIF713:       ENDIF
767: !714: !
768: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.715: ! Skip image INTIMAGE+2 - no non-zero gradients on other images and no energy contributions.
769: !716: !
770: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1717: !     IF ((D2.LT.INTCONSTRAINREPCUT).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
771:  
772: !     IF ((NREPJ(J2).EQ.2024).OR.(NREPI(J2).EQ.2024)) WRITE(*, '(A,2I6,3G20.10,I6)') & 
773: ! &    'j,i,D2,NREPCUT,INTCONST,DUMMY,J1=',NREPJ(J2),NREPI(J2),D2,NREPCUT(J2),INTCONST,J1 
774:  
775:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1718:       IF ((D2.LT.NREPCUT(J2)).AND.(J1.LT.INTIMAGE+2)) THEN ! terms for image J1 - non-zero derivatives only for J1
776: !        D12=DSQ2**6719: !        D12=DSQ2**6
777:          D12=DSQ2720:          D12=DSQ2
778: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)721: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
779: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)722: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D2-13.0D0*NREPCUT(J2))/INTCONST)
780:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)723:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D2-3.0D0*NREPCUT(J2))/INTCONST)
781: !        IF ((NREPJ(J2).EQ.2024).OR.(NREPI(J2).EQ.2024)) WRITE(*, '(A,2I6,6G20.10)') & 
782: ! &        'j,i,INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',NREPJ(J2),NREPI(J2),INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY  
783: !        IF ((NREPJ(J2).EQ.83).AND.(NREPI(J2).EQ.357)) WRITE(*, '(A,6G20.10)') &724: !        IF ((NREPJ(J2).EQ.83).AND.(NREPI(J2).EQ.357)) WRITE(*, '(A,6G20.10)') &
784: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 725: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 
785: !        IF ((NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.357)) WRITE(*, '(A,6G20.10)') &726: !        IF ((NREPI(J2).EQ.83).AND.(NREPJ(J2).EQ.357)) WRITE(*, '(A,6G20.10)') &
786: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 727: ! &               'INTCONSTRAINTREP,D12,D2,NREPCUT,INTCONST,DUMMY=',INTCONSTRAINTREP,D12,D2,NREPCUT(J2),INTCONST ,DUMMY 
787:          EEE(J1)=EEE(J1)+DUMMY728:          EEE(J1)=EEE(J1)+DUMMY
788:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN729:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
789:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &730:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R1 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
790:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY731:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
791:          ENDIF732:          ENDIF
792:          IF (DUMMY.GT.EMAX) THEN 
793:             IMAX=J1 
794:             JMAX=J2 
795:             EMAX=DUMMY 
796:          ENDIF 
797:         IF (DUMMY.GT.EEMAX(J1)) THEN 
798:             JJMAX(J1)=J2 
799:             EEMAX(J1)=DUMMY 
800:          ENDIF 
801:          REPE(J1)=REPE(J1)+DUMMY733:          REPE(J1)=REPE(J1)+DUMMY
802:          EREP=EREP+DUMMY734:          EREP=EREP+DUMMY
803: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)735: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
804:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)736:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D2*D12)-1.0D0/INTCONST)
805:          REPGRAD(1:3)=DUMMY*G2(1:3)737:          REPGRAD(1:3)=DUMMY*G2(1:3)
806: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &738: !        WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
807: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)739: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
808:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)740:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
809:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)741:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
810:       ENDIF742:       ENDIF
817: !        D12=DSQ1**6749: !        D12=DSQ1**6
818:          D12=DSQ1750:          D12=DSQ1
819: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*INTCONSTRAINREPCUT)/INTCONST)751: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*INTCONSTRAINREPCUT)/INTCONST)
820: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*NREPCUT(J2))/INTCONST)752: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*D1-13.0D0*NREPCUT(J2))/INTCONST)
821:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D1-3.0D0*NREPCUT(J2))/INTCONST)753:          DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*D1-3.0D0*NREPCUT(J2))/INTCONST)
822:          EEE(J1-1)=EEE(J1-1)+DUMMY754:          EEE(J1-1)=EEE(J1-1)+DUMMY
823:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN755:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
824:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R2 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &756:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R2 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
825:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY757:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
826:          ENDIF758:          ENDIF
827:          IF (DUMMY.GT.EMAX) THEN 
828:             IMAX=J1 
829:             JMAX=J2 
830:             EMAX=DUMMY 
831:          ENDIF 
832:          IF (DUMMY.GT.EEMAX(J1-1)) THEN 
833:             JJMAX(J1-1)=J2 
834:             EEMAX(J1-1)=DUMMY 
835:          ENDIF 
836:          REPE(J1-1)=REPE(J1-1)+DUMMY759:          REPE(J1-1)=REPE(J1-1)+DUMMY
837:          EREP=EREP+DUMMY760:          EREP=EREP+DUMMY
838: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)761: !        DUMMY=-12.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)
839:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)762:          DUMMY=-2.0D0*INTCONSTRAINTREP*(1.0D0/(D1*D12)-1.0D0/INTCONST)
840:          REPGRAD(1:3)=DUMMY*G1(1:3)763:          REPGRAD(1:3)=DUMMY*G1(1:3)
841: !        WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &764: !        WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,REPI,REPJ,REPGRAD=',J1,J2,NREPI(J2),NREPJ(J2), &
842: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)765: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
843:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)766:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
844:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)767:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
845:       ENDIF768:       ENDIF
852: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)775: !        DUMMY=INTCONSTRAINTREP*(1.0D0/D12+(12.0D0*DINT-13.0D0*NREPCUT(J2))/INTCONST)
853:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)776:          DUMMY=INTMINFAC*INTCONSTRAINTREP*(1.0D0/D12+(2.0D0*DINT-3.0D0*NREPCUT(J2))/INTCONST)
854:          EREP=EREP+DUMMY777:          EREP=EREP+DUMMY
855: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &778: !        IF (DUMMY.GT.1.0D7) PRINT '(A,3I6,3G20.10)','J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY=', &
856: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY779: ! &                                                   J2,NREPI(J2),NREPJ(J2),DINT,NREPCUT(J2),DUMMY
857: !        ENDIF780: !        ENDIF
858:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN781:          IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
859:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &782:             WRITE(*, '(A,2I6,2L5,G20.10)') 'R3 NREPI,NREPJ,INTFROZEN(NREPI),INTFROZEN(NREPJ),DUMMY=', &
860:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY783:   &                                     NREPI(J2),NREPJ(J2),INTFROZEN(NREPI(J2)),INTFROZEN(NREPJ(J2)),DUMMY
861:          ENDIF784:          ENDIF
862:          IF (DUMMY.GT.EMAX) THEN 
863:             IMAX=J1 
864:             JMAX=J2 
865:             EMAX=DUMMY 
866:          ENDIF 
867:          IF (DUMMY.GT.EEMAX(J1)) THEN 
868:             JJMAX(J1)=J2 
869:             EEMAX(J1)=DUMMY 
870:          ENDIF 
871:          IF (DUMMY.GT.EEMAX(J1-1)) THEN 
872:             JJMAX(J1-1)=J2 
873:             EEMAX(J1-1)=DUMMY 
874:          ENDIF 
875:          IF (J1.EQ.2) THEN785:          IF (J1.EQ.2) THEN
876:             EEE(J1)=EEE(J1)+DUMMY786:             EEE(J1)=EEE(J1)+DUMMY
877:             REPEINT(J1)=REPEINT(J1)+DUMMY787:             REPEINT(J1)=REPEINT(J1)+DUMMY
878:             NREPINT(J1)=NREPINT(J1)+1788:             NREPINT(J1)=NREPINT(J1)+1
879:          ELSE IF (J1.LT.INTIMAGE+2) THEN789:          ELSE IF (J1.LT.INTIMAGE+2) THEN
880:             EEE(J1)=EEE(J1)+DUMMY/2.0D0790:             EEE(J1)=EEE(J1)+DUMMY/2.0D0
881:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0791:             EEE(J1-1)=EEE(J1-1)+DUMMY/2.0D0
882:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0792:             REPEINT(J1)=REPEINT(J1)+DUMMY/2.0D0
883:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0793:             REPEINT(J1-1)=REPEINT(J1-1)+DUMMY/2.0D0
884:             NREPINT(J1)=NREPINT(J1)+1794:             NREPINT(J1)=NREPINT(J1)+1
894: !        WRITE(*, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &804: !        WRITE(*, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
895: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)805: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
896:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)806:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
897:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)807:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
898:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)808:          REPGRAD(1:3)=INTMINFAC*DUMMY*G2INT(1:3)
899: !        WRITE(*, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &809: !        WRITE(*, '(A,4I6,2G15.5)') 'in1 J1,J2,REPI,REPJ,REPGRAD,NREPCUT=',J1,J2,NREPI(J2),NREPJ(J2), &
900: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)810: ! &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2),NREPCUT(J2)
901:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)811:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
902:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)812:          GGG(NJ2+1:NJ2+3)=GGG(NJ2+1:NJ2+3)-REPGRAD(1:3)
903:       ENDIF813:       ENDIF
 814: !     WRITE(*,'(A,2I8,6G20.10)') 'congrad2> C J1,J2,GGG(1:6)=',J1,J2,GGG(1:6)
904:    ENDDO815:    ENDDO
905: ENDDO816: ENDDO
906: ! DO J2=1,INTIMAGE+2 
907: !    IF (JJMAX(J2).GT.0) THEN 
908: !       WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest repulsive  contribution for image ',J2,' pair index ', & 
909: !  &                                JJMAX(J2),' atoms ', & 
910: !  &                                NREPI(JJMAX(J2)),NREPJ(JJMAX(J2)) 
911: !    ENDIF 
912: ! ENDDO 
913: IF (JMAX.GT.0) THEN 
914:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest repulsive  contribution for any image in image ',IMAX, & 
915:  &  ' pair index ', & 
916:  &                                JMAX,' atoms ',NREPI(JMAX),NREPJ(JMAX) 
917: ENDIF 
918: !817: !
919: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise818: ! Spring energy. Set EEE(J1) and ESPRING dividing up the pairwise
920: ! energy terms between images except for the end points.819: ! energy terms between images except for the end points.
921: !820: !
922: ESPRING=0.0D0821: ESPRING=0.0D0
923: EMAX=0.0D0 
924: IMAX=0 
925: IF (KINT.NE.0.0D0) THEN822: IF (KINT.NE.0.0D0) THEN
926:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1823:    DO J1=1,INTIMAGE+1 ! sum over edges from J1 to J1+1
927:       NI1=(3*NATOMS)*(J1-1)824:       NI1=(3*NATOMS)*(J1-1)
928:       NI2=(3*NATOMS)*J1825:       NI2=(3*NATOMS)*J1
929: !826: !
930: !  Edge between J1 and J1+1827: !  Edge between J1 and J1+1
931: !828: !
932:       DPLUS=0.0D0829:       DPLUS=0.0D0
933:       DO J2=1,NATOMS830:       DO J2=1,NATOMS
934:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 831:          IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
935:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &832:             DPLUS=DPLUS+(XYZ(NI1+3*(J2-1)+1)-XYZ(NI2+3*(J2-1)+1))**2 &
936:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &833:   &                    +(XYZ(NI1+3*(J2-1)+2)-XYZ(NI2+3*(J2-1)+2))**2 &
937:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2834:   &                    +(XYZ(NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+3))**2
938:          ENDIF835:          ENDIF
939: !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS836: !        WRITE(*,'(A,2I8,G20.10)') 'J1,J2,DPLUS: ',J1,J2,DPLUS
940:       ENDDO837:       ENDDO
941:       DPLUS=SQRT(DPLUS)838:       DPLUS=SQRT(DPLUS)
942:       IF (DPLUS.GT.IMSEPMAX) THEN839:       IF (DPLUS.GT.IMSEPMAX) THEN
943:          DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2840:          DUMMY=KINT*0.5D0*(DPLUS-IMSEPMAX)**2
944:          IF (DUMMY.GT.EMAX) THEN 
945:             IMAX=J1 
946:             EMAX=DUMMY 
947:          ENDIF 
948: !        DUMMY=KINT*0.5D0*DPLUS**2841: !        DUMMY=KINT*0.5D0*DPLUS**2
949:          ESPRING=ESPRING+DUMMY842:          ESPRING=ESPRING+DUMMY
950:          IF (DUMMY.GT.EEMAX(J1+1)) THEN 
951:             EEMAX(J1+1)=DUMMY 
952:          ENDIF 
953:  
954: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING843: !        WRITE(*,'(A,4G20.10)') 'DPLUS,IMSEPMAX,DUMMY,ESPRING=',DPLUS,IMSEPMAX,DUMMY,ESPRING
955:          DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS844:          DUMMY=KINT*(DPLUS-IMSEPMAX)/DPLUS
956: !        DUMMY=KINT845: !        DUMMY=KINT
957:          DO J2=1,NATOMS846:          DO J2=1,NATOMS
958:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 847:             IF ((.NOT.INTSPRINGACTIVET).OR.ATOMACTIVE(J2)) THEN 
959:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))848:                SPGRAD(1:3)=DUMMY*(XYZ(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)-XYZ(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3))
960:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)849:                GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)=GGG(NI1+3*(J2-1)+1:NI1+3*(J2-1)+3)+SPGRAD(1:3)
961:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3)850:                GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)=GGG(NI2+3*(J2-1)+1:NI2+3*(J2-1)+3)-SPGRAD(1:3)
962:             ENDIF851:             ENDIF
963:          ENDDO852:          ENDDO
964:       ENDIF853:       ENDIF
965:    ENDDO854:    ENDDO
966: ENDIF855: ENDIF
967: WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad> Highest spring  contribution for any image in image ',IMAX 
968:          IF (PRINTE) THEN856:          IF (PRINTE) THEN
969:             WRITE(*, '(A,G20.10)') 'ESPRING=',ESPRING857:             WRITE(*, '(A,G20.10)') 'ESPRING=',ESPRING
970:          ENDIF858:          ENDIF
971: !859: !
972: ! Set gradients on frozen atoms to zero.860: ! Set gradients on frozen atoms to zero.
973: !861: !
974: IF (FREEZE) THEN862: IF (FREEZE) THEN
975:    DO J1=2,INTIMAGE+1  863:    DO J1=2,INTIMAGE+1  
976:       DO J2=1,NATOMS864:       DO J2=1,NATOMS
977:          IF (FROZEN(J2)) THEN865:          IF (FROZEN(J2)) THEN


r33407/intlbfgs.f90 2017-10-23 10:30:21.257310020 +0100 r33406/intlbfgs.f90 2017-10-23 10:30:22.945332350 +0100
 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, & 23:      & INTCONSTRAINREPCUT, REPCON, INTCONSTRAINTREP, INTREPSEP, NREPI, NREPJ, &
 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, & 24:      & CONDISTREFLOCAL, INTCONFRAC, CONACTIVE, REPI, &
 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, & 25:      & REPJ, NREPMAX, ATOMACTIVE, NCONSTRAINTON, CONION, CONJON, CONDISTREFLOCALON, CONDISTREFON, &
 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, & 26:      & NREPCUT, REPCUT, CHECKCONINT, INTCONSTEPS, INTRELSTEPS, MAXCONE, COLDFUSIONLIMIT, &
 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, & 27:      & INTSTEPS1, DUMPINTXYZ, DUMPINTXYZFREQ, DUMPINTEOS, DUMPINTEOSFREQ, &
 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, & 28:      & IMSEPMIN, IMSEPMAX, MAXINTIMAGE, INTFREEZET, INTFREEZETOL, FREEZE, &
 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, & 29:      & INTFROZEN, CHECKREPINTERVAL, NNREPULSIVE, INTFREEZEMIN, INTIMAGECHECK, &
 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, & 30:      & CONCUT, CONCUTLOCAL, KINT, REPIFIX, REPJFIX, NREPULSIVEFIX, &
 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, & 31:      & NCONSTRAINTFIX, CONIFIX, CONJFIX, QCIPERMCHECK, QCIPERMCHECKINT, BULKT, TWOD, RIGIDBODY, &
 32:      & QCIADDREP, QCIXYZ, WHOLEDNEB, QCIIMAGE, FROZEN, QCIRESTART, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, & 32:      & QCIADDREP, QCIXYZ, WHOLEDNEB, QCIIMAGE, FROZEN, QCIRESTART, NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS, &
 33:      & PERMDIST, LOCALPERMCUT, QCILPERMDIST, QCIPDINT, QCIPERMCUT, QCIAMBERT, BONDS 33:      & PERMDIST, LOCALPERMCUT, QCILPERMDIST, QCIPDINT, QCIPERMCUT
 34: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3 34: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3
 35: USE MODCHARMM, ONLY : CHRMMT 35: USE MODCHARMM, ONLY : CHRMMT
 36: USE CHIRALITY 
 37:  36: 
 38: IMPLICIT NONE  37: IMPLICIT NONE 
 39:  38: 
 40: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points 39: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points
 41: INTEGER D, U 40: INTEGER D, U
 42: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE, NEIGHBOUR_COORDS(12), CENTRE_COORDS(3) 41: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE
 43: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ 42: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ
 44: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF   43: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF  
 45: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS), STARTGROUP(NPERMGROUP), ENDGROUP(NPERMGROUP) 44: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS)
 46: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS), IDENTITY 45: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS)
 47: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 46: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 48:  47: 
 49: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY 48: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY
 50: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF 49: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF
 51: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2 50: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2
 52: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND 51: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE
 53: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS)  52: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS)
 54: LOGICAL CHIRALSR, CHIRALSRP  
 55: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS) 53: INTEGER NTRIES(NATOMS), NITERDONE, EXITSTATUS, DLIST(NATOMS)
 56: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, & 54: DOUBLE PRECISION :: DDOT,STPMIN, ETOTALTMP, RMSTMP, USEFRAC, STIME, FTIME, &
 57:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), & 55:   &                 ETOTAL, LASTGOODE, RMS, STEPTOT, LINTCONSTRAINTTOL, LXYZ(2*3*NATOMS), &
 58:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS) 56:   &                 BESTWORST, WORST, COORDSA(3*NATOMS), COORDSB(3*NATOMS), COORDSC(3*NATOMS)
 59: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA 57: DOUBLE PRECISION, DIMENSION(MUPDATE)     :: RHO1,ALPHA
 60: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS) 58: DOUBLE PRECISION :: EOLD, DMOVED(NATOMS)
 61: LOGICAL SWITCHED 59: LOGICAL SWITCHED
 62: DOUBLE PRECISION, POINTER :: X(:), G(:) 60: DOUBLE PRECISION, POINTER :: X(:), G(:)
 63: ! 61: !
 64: ! efk: for freezenodes 62: ! efk: for freezenodes
 68: ! 66: !
 69: ! Dimensions involving INTIMAGE 67: ! Dimensions involving INTIMAGE
 70: ! 68: !
 71: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), & 69: DOUBLE PRECISION, ALLOCATABLE :: TRUEEE(:), &
 72:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), & 70:   &              EEETMP(:), MYGTMP(:), EEE(:), STEPIMAGE(:), &
 73:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:) 71:   &              GTMP(:), DIAG(:), STP(:), SEARCHSTEP(:,:), GDIF(:,:), GLAST(:), XSAVE(:)
 74: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:) 72: DOUBLE PRECISION, ALLOCATABLE, TARGET :: XYZ(:), GGG(:), DPTMP(:), D2TMP(:,:)
 75: ! saved interpolation 73: ! saved interpolation
 76: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE 74: INTEGER BESTINTIMAGE, NSTEPS, NITERUSE
 77: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:) 75: LOGICAL, ALLOCATABLE :: CHECKG(:), IMGFREEZE(:)
 78: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP), CHIRALACTIVE(NATOMS) 76: LOGICAL READIMAGET, GROUPACTIVE(NPERMGROUP)
 79: INTEGER LUNIT, GETUNIT 77: INTEGER LUNIT, GETUNIT
 80: CHARACTER(LEN=2) SDUMMY 78: CHARACTER(LEN=2) SDUMMY
 81: INTEGER JMAXEEE,JMAXRMS 79: INTEGER JMAXEEE,JMAXRMS
 82: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT 80: DOUBLE PRECISION MAXEEE,MAXRMS,MINEEE,SAVELOCALPERMCUT
 83:  81: 
 84: WHOLEDNEB=.FALSE. 82: WHOLEDNEB=.FALSE.
 85: READIMAGET=.FALSE. 83: READIMAGET=.FALSE.
 86: REMOVEIMAGE=.FALSE. 84: REMOVEIMAGE=.FALSE.
 87: IF (QCIRESTART) READIMAGET=.TRUE. 85: IF (QCIRESTART) READIMAGET=.TRUE.
 88: IF (READIMAGET) THEN 86: IF (READIMAGET) THEN
120: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',MUPDATE,' maximum step size=',MAXBFGS118: WRITE(*,'(A,I6,A,G20.10)') ' intlbfgs> Updates: ',MUPDATE,' maximum step size=',MAXBFGS
121: ADDATOM=.FALSE.119: ADDATOM=.FALSE.
122: NFAIL=0120: NFAIL=0
123: IMGFREEZE(1:INTIMAGE)=.FALSE.121: IMGFREEZE(1:INTIMAGE)=.FALSE.
124: D=(3*NATOMS)*INTIMAGE122: D=(3*NATOMS)*INTIMAGE
125: U=MUPDATE123: U=MUPDATE
126: NITERDONE=1124: NITERDONE=1
127: NITERUSE=1125: NITERUSE=1
128: NQDONE=0126: NQDONE=0
129: 127: 
130: ! 
131: ! Read AMBER topology and count bonded neighbours for each atom. 
132: ! Assume that we won't have more than 6 bonds! 
133: ! 
134: IF (QCIAMBERT) THEN 
135:    CALL TOPOLOGY_READER(NBOND) 
136:    NBONDED(1:NATOMS)=0 
137:    BONDEDLIST(1:NATOMS,1:6)=0 
138:    DO J2=1,NBOND     
139:       NBONDED(BONDS(J2,1))=NBONDED(BONDS(J2,1))+1 
140:       IF (NBONDED(BONDS(J2,1)).GT.6) THEN 
141:          WRITE(*,'(A,I6)') 'intlbfgs> ERROR *** atom ',BONDS(J2,1),' has more than 6 bonds' 
142:          STOP 
143:       ENDIF 
144:       BONDEDLIST(BONDS(J2,1),NBONDED(BONDS(J2,1)))=BONDS(J2,2) 
145:       NBONDED(BONDS(J2,2))=NBONDED(BONDS(J2,2))+1 
146:       IF (NBONDED(BONDS(J2,2)).GT.6) THEN 
147:          WRITE(*,'(A,I6)') 'intlbfgs> ERROR *** atom ',BONDS(J2,2),' has more than 6 bonds' 
148:          STOP 
149:       ENDIF 
150:       BONDEDLIST(BONDS(J2,2),NBONDED(BONDS(J2,2)))=BONDS(J2,1) 
151:    ENDDO 
152: !  DO J2=1,NATOMS 
153: !     WRITE(*,'(A,I6,A,I6,A,6I6)') 'intlbfgs> atom ',J2,' bonds ',NBONDED(J2),' to atoms ',BONDEDLIST(J2,1:NBONDED(J2)) 
154: !  ENDDO 
155: ENDIF 
156:  
157: IF ( D<=0 ) THEN128: IF ( D<=0 ) THEN
158:    WRITE(*,*) 'd is not positive, d=',d129:    WRITE(*,*) 'd is not positive, d=',d
159:    STOP130:    STOP
160: ENDIF131: ENDIF
161: IF ( U<=0 ) THEN132: IF ( U<=0 ) THEN
162:    WRITE(*,*) 'u is not positive, u=',u133:    WRITE(*,*) 'u is not positive, u=',u
163:    STOP134:    STOP
164: ENDIF135: ENDIF
165: IF (INTSTEPS1 < 0) THEN136: IF (INTSTEPS1 < 0) THEN
166:    WRITE(*,'(1x,a)') 'Maximal number of iterations is less than zero! Stop.'137:    WRITE(*,'(1x,a)') 'Maximal number of iterations is less than zero! Stop.'
172: !143: !
173: IF (.NOT.ALLOCATED(CONI)) THEN 144: IF (.NOT.ALLOCATED(CONI)) THEN 
174:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX))145:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX))
175:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))146:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
176: ENDIF147: ENDIF
177: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))148: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
178: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))149: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
179: !150: !
180: ! Initialise XYZ151: ! Initialise XYZ
181: !152: !
182: IF (READIMAGET) THEN  ! Note that this will ignore the coordinates in start and finish153: IF (READIMAGET) THEN
183:    LUNIT=GETUNIT()154:    LUNIT=GETUNIT()
184:    OPEN(UNIT=LUNIT,FILE='int.xyz',STATUS='OLD')155:    OPEN(UNIT=LUNIT,FILE='int.xyz',STATUS='OLD')
185:    DO J2=1,INTIMAGE+2156:    DO J2=1,INTIMAGE+2
186:       READ(LUNIT,*) NDUMMY157:       READ(LUNIT,*) NDUMMY
187:       READ(LUNIT,*) 158:       READ(LUNIT,*) 
188:       DO J1=1,NATOMS159:       DO J1=1,NATOMS
189:          READ(LUNIT,*) SDUMMY,XYZ(3*NATOMS*(J2-1)+3*(J1-1)+1),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+2),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+3)160:          READ(LUNIT,*) SDUMMY,XYZ(3*NATOMS*(J2-1)+3*(J1-1)+1),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+2),XYZ(3*NATOMS*(J2-1)+3*(J1-1)+3)
190:       ENDDO161:       ENDDO
191:    ENDDO162:    ENDDO
192:    CLOSE(LUNIT)163:    CLOSE(LUNIT)
193: ! 
194: ! Don't overwrite start and finish - we should have aligned these by now 
195: ! 
196:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS)) 
197:    XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS)) 
198: ELSE164: ELSE
199:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))165:    XYZ(1:(3*NATOMS))=QSTART(1:(3*NATOMS))
200:    XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS))166:    XYZ((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=QFINISH(1:(3*NATOMS))
201:    DO J1=1,INTIMAGE+2167:    DO J1=1,INTIMAGE+2
202:       XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)168:       XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)
203:    ENDDO169:    ENDDO
204: ENDIF170: ENDIF
205: 171: 
206: NQCIFREEZE=0172: NQCIFREEZE=0
207: ! IF (FREEZE) THEN173: ! IF (FREEZE) THEN
281:    DO J1=1,INTIMAGE+2247:    DO J1=1,INTIMAGE+2
282:       XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)248:       XYZ((J1-1)*(3*NATOMS)+1:J1*(3*NATOMS))=((INTIMAGE+2-J1)*QSTART(1:(3*NATOMS))+(J1-1)*QFINISH(1:(3*NATOMS)))/(INTIMAGE+1)
283:    ENDDO249:    ENDDO
284:    GOTO 678250:    GOTO 678
285: ENDIF251: ENDIF
286: 252: 
287: NACTIVE=0253: NACTIVE=0
288: TURNONORDER(1:NATOMS)=0254: TURNONORDER(1:NATOMS)=0
289: ATOMACTIVE(1:NATOMS)=.FALSE.255: ATOMACTIVE(1:NATOMS)=.FALSE.
290: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT**6 ! also needed for congrad.f90 potential256: REPCON=-INTCONSTRAINTREP/INTCONSTRAINREPCUT**6 ! also needed for congrad.f90 potential
291: IF (QCIRESTART) THEN257: IF (READIMAGET) THEN
292:    LUNIT=GETUNIT()258:    LUNIT=GETUNIT()
293:    OPEN(UNIT=LUNIT,FILE='QCIdump',STATUS='OLD')259:    OPEN(UNIT=LUNIT,FILE='QCIdump',STATUS='OLD')
294: 260: 
295:    READ(LUNIT,*) NACTIVE261:    READ(LUNIT,*) NACTIVE
296:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NACTIVE,' active atoms'262:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NACTIVE,' active atoms'
297:    WRITE(*,'(A)') ' intlbfgs> reading turnonorder'263:    WRITE(*,'(A)') ' intlbfgs> reading turnonorder'
298:    READ(LUNIT,*) TURNONORDER(1:NACTIVE)264:    READ(LUNIT,*) TURNONORDER(1:NACTIVE)
299: !  WRITE(*,'(12I8)') TURNONORDER(1:NACTIVE)265:    WRITE(*,'(12I8)') TURNONORDER(1:NACTIVE)
300:    WRITE(*,'(A)') ' intlbfgs> reading atomactive'266:    WRITE(*,'(A)') ' intlbfgs> reading atomactive'
301:    READ(LUNIT,*) ATOMACTIVE(1:NATOMS)267:    READ(LUNIT,*) ATOMACTIVE(1:NATOMS)
302: !  WRITE(*,'(12L5)') ATOMACTIVE(1:NATOMS) 268:    WRITE(*,'(12L5)') ATOMACTIVE(1:NATOMS) 
303:    READ(LUNIT,*) NCONSTRAINT269:    READ(LUNIT,*) NCONSTRAINT
304:    WRITE(*,'(A)') ' intlbfgs> reading conactive'270:    WRITE(*,'(A)') ' intlbfgs> reading conactive'
305:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NCONSTRAINT,' constraints'271:    WRITE(*,'(A,I10,A)') ' intlbfgs> restart has ',NCONSTRAINT,' constraints'
306:    ALLOCATE(CONACTIVE(NCONSTRAINT))272:    ALLOCATE(CONACTIVE(NCONSTRAINT))
307:    READ(LUNIT,*) CONACTIVE(1:NCONSTRAINT)273:    READ(LUNIT,*) CONACTIVE(1:NCONSTRAINT)
308: !  WRITE(*,'(12L5)') CONACTIVE(1:NCONSTRAINT) 274:    WRITE(*,'(12L5)') CONACTIVE(1:NCONSTRAINT) 
309: 275: 
310:    ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))276:    ALLOCATE(CONDISTREFLOCAL(NCONSTRAINT))
311:    ALLOCATE(CONCUTLOCAL(NCONSTRAINT))277:    ALLOCATE(CONCUTLOCAL(NCONSTRAINT))
312:    CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)278:    CONDISTREFLOCAL(1:NCONSTRAINT)=CONDISTREF(1:NCONSTRAINT)
313:    CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)279:    CONCUTLOCAL(1:NCONSTRAINT)=CONCUT(1:NCONSTRAINT)
314: !  WRITE(*,'(A)') ' intlbfgs> CONCUT:'280:    WRITE(*,'(A)') ' intlbfgs> CONCUT:'
315: !  WRITE(*,'(6G20.10)') CONCUT(1:NCONSTRAINT)281:    WRITE(*,'(6G20.10)') CONCUT(1:NCONSTRAINT)
316: !  WRITE(*,'(A)') ' intlbfgs> CONDISTREF:'282:    WRITE(*,'(A)') ' intlbfgs> CONDISTREF:'
317: !  WRITE(*,'(6G20.10)') CONDISTREF(1:NCONSTRAINT)283:    WRITE(*,'(6G20.10)') CONDISTREF(1:NCONSTRAINT)
318: !  WRITE(*,'(A)') ' intlbfgs> CONI:'284:    WRITE(*,'(A)') ' intlbfgs> CONI:'
319: !  WRITE(*,'(12I8)') CONI(1:NCONSTRAINT)285:    WRITE(*,'(12I8)') CONI(1:NCONSTRAINT)
320: !  WRITE(*,'(A)') ' intlbfgs> CONJ:'286:    WRITE(*,'(A)') ' intlbfgs> CONJ:'
321: !  WRITE(*,'(12I8)') CONJ(1:NCONSTRAINT)287:    WRITE(*,'(12I8)') CONJ(1:NCONSTRAINT)
322: 288: 
323:    READ(LUNIT,*) NREPULSIVE,NNREPULSIVE,NREPMAX289:    READ(LUNIT,*) NREPULSIVE,NNREPULSIVE,NREPMAX
324:    USEFRAC=1.0D0290:    USEFRAC=1.0D0
325: !  WRITE(*,'(A,3I10,G20.10)') 'intlbfgs> NREPULSIVE,NNREPULSIVE,NREPMAX=',NREPULSIVE,NNREPULSIVE,NREPMAX291:    WRITE(*,'(A,3I10,G20.10)') 'intlbfgs> NREPULSIVE,NNREPULSIVE,NREPMAX=',NREPULSIVE,NNREPULSIVE,NREPMAX
326:    IF (ALLOCATED(REPI)) DEALLOCATE(REPI)292:    IF (ALLOCATED(REPI)) DEALLOCATE(REPI)
327:    IF (ALLOCATED(REPJ)) DEALLOCATE(REPJ)293:    IF (ALLOCATED(REPJ)) DEALLOCATE(REPJ)
328:    IF (ALLOCATED(NREPI)) DEALLOCATE(NREPI)294:    IF (ALLOCATED(NREPI)) DEALLOCATE(NREPI)
329:    IF (ALLOCATED(NREPJ)) DEALLOCATE(NREPJ)295:    IF (ALLOCATED(NREPJ)) DEALLOCATE(NREPJ)
330:    IF (ALLOCATED(REPCUT)) DEALLOCATE(REPCUT)296:    IF (ALLOCATED(REPCUT)) DEALLOCATE(REPCUT)
331:    IF (ALLOCATED(NREPCUT)) DEALLOCATE(NREPCUT)297:    IF (ALLOCATED(NREPCUT)) DEALLOCATE(NREPCUT)
332:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))298:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
333:    READ(LUNIT,*) REPI(1:NREPULSIVE)299:    READ(LUNIT,*) REPI(1:NREPULSIVE)
334:    WRITE(*,'(A)') ' intlbfgs> read REPI:'300:    WRITE(*,'(A)') ' intlbfgs> read REPI:'
335: !  WRITE(*,'(12I8)') REPI(1:NREPULSIVE)301:    WRITE(*,'(12I8)') REPI(1:NREPULSIVE)
336:    READ(LUNIT,*) REPJ(1:NREPULSIVE)302:    READ(LUNIT,*) REPJ(1:NREPULSIVE)
337:    WRITE(*,'(A)') ' intlbfgs> read REPJ:'303:    WRITE(*,'(A)') ' intlbfgs> read REPJ:'
338: !  WRITE(*,'(12I8)') REPJ(1:NREPULSIVE)304:    WRITE(*,'(12I8)') REPJ(1:NREPULSIVE)
339:    READ(LUNIT,*) NREPI(1:NNREPULSIVE)305:    READ(LUNIT,*) NREPI(1:NNREPULSIVE)
340:    WRITE(*,'(A)') ' intlbfgs> read NREPI:'306:    WRITE(*,'(A)') ' intlbfgs> read NREPI:'
341: !  WRITE(*,'(12I8)') NREPI(1:NNREPULSIVE)307:    WRITE(*,'(12I8)') NREPI(1:NNREPULSIVE)
342:    READ(LUNIT,*) NREPJ(1:NNREPULSIVE)308:    READ(LUNIT,*) NREPJ(1:NNREPULSIVE)
343:    WRITE(*,'(A)') ' intlbfgs> read NREPJ:'309:    WRITE(*,'(A)') ' intlbfgs> read NREPJ:'
344: !  WRITE(*,'(12I8)') NREPJ(1:NNREPULSIVE)310:    WRITE(*,'(12I8)') NREPJ(1:NNREPULSIVE)
345:    READ(LUNIT,*) REPCUT(1:NREPULSIVE)311:    READ(LUNIT,*) REPCUT(1:NREPULSIVE)
346:    WRITE(*,'(A)') ' intlbfgs> read REPCUT:'312:    WRITE(*,'(A)') ' intlbfgs> read REPCUT:'
347: !  WRITE(*,'(6G20.10)') REPCUT(1:NREPULSIVE)313:    WRITE(*,'(6G20.10)') REPCUT(1:NREPULSIVE)
348:    READ(LUNIT,*) NREPCUT(1:NNREPULSIVE)314:    READ(LUNIT,*) NREPCUT(1:NNREPULSIVE)
349:    WRITE(*,'(A)') ' intlbfgs> read NREPCUT:'315:    WRITE(*,'(A)') ' intlbfgs> read NREPCUT:'
350: !  WRITE(*,'(6G20.10)') NREPCUT(1:NNREPULSIVE)316:    WRITE(*,'(6G20.10)') NREPCUT(1:NNREPULSIVE)
351:    READ(LUNIT,*) INTFROZEN(1:NATOMS)317:    READ(LUNIT,*) INTFROZEN(1:NATOMS)
352:    WRITE(*,'(A)') ' intlbfgs> read INTFROZEN'318:    WRITE(*,'(A)') ' intlbfgs> read INTFROZEN'
353: !  WRITE(*,'(12L5)') INTFROZEN(1:NATOMS)319:    WRITE(*,'(12L5)') INTFROZEN(1:NATOMS)
354:    CLOSE(LUNIT)320:    CLOSE(LUNIT)
355: 321: 
356:    GLAST(1:D)=G(1:D)322:    GLAST(1:D)=G(1:D)
357:    XSAVE(1:D)=X(1:D)323:    XSAVE(1:D)=X(1:D)
358:    GOTO 986324:    GOTO 986
359: ENDIF325: ENDIF
360: IF (INTFREEZET) THEN326: IF (INTFREEZET) THEN
361:    DO J1=1,NATOMS327:    DO J1=1,NATOMS
362:       IF (INTFROZEN(J1)) THEN328:       IF (INTFROZEN(J1)) THEN
363: ! 329: ! 
575: 541: 
576: !542: !
577: ! In this block PERMGROUP(NDUMMY+J2-1) counts through the atom indices specified in perm.allow,543: ! In this block PERMGROUP(NDUMMY+J2-1) counts through the atom indices specified in perm.allow,
578: ! one group at a time.544: ! one group at a time.
579: !545: !
580: IF (PERMDIST) THEN546: IF (PERMDIST) THEN
581:    PERMUTABLE(1:NATOMS)=.FALSE.547:    PERMUTABLE(1:NATOMS)=.FALSE.
582:    INGROUP(1:NATOMS)=0548:    INGROUP(1:NATOMS)=0
583:    NDUMMY=1549:    NDUMMY=1
584:    DO J1=1,NPERMGROUP550:    DO J1=1,NPERMGROUP
585:       STARTGROUP(J1)=NDUMMY 
586:       GROUPACTIVE(J1)=.FALSE.551:       GROUPACTIVE(J1)=.FALSE.
587:       DO J2=1,NPERMSIZE(J1)552:       DO J2=1,NPERMSIZE(J1)
588:          PERMUTABLE(PERMGROUP(NDUMMY+J2-1))=.TRUE.553:          PERMUTABLE(PERMGROUP(NDUMMY+J2-1))=.TRUE.
589:          INGROUP(PERMGROUP(NDUMMY+J2-1))=J1554:          INGROUP(PERMGROUP(NDUMMY+J2-1))=J1
590:       ENDDO555:       ENDDO
591:       NDUMMY=NDUMMY+NPERMSIZE(J1)556:       NDUMMY=NDUMMY+NPERMSIZE(J1)
592:       ENDGROUP(J1)=NDUMMY-1 
593:    ENDDO557:    ENDDO
594: ENDIF558: ENDIF
595: 559: 
596: 567 CONTINUE560: 567 CONTINUE
597: 561: 
598: DO ! Main do loop with counter NITERDONE, initially set to one562: DO ! Main do loop with counter NITERDONE, initially set to one
599: 563: 
600: !564: !
601: !  Check permutational alignments. Maintain a list of the permutable groups where all565: !  Check permutational alignments. Maintain a list of the permutable groups where all
602: !  members are active. See if we have any new complete groups. MUST update NDUMMY566: !  members are active. See if we have any new complete groups. MUST update NDUMMY
603: !  counter to step through permutable atom list.567: !  counter to step through permutable atom list.
604: !568: !
605: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN569: IF (QCILPERMDIST.AND.(MOD(NITERDONE-1,QCIPDINT).EQ.0)) THEN
606:  
607:    CHIRALACTIVE(1:NATOMS)=.FALSE. 
608:    chicheck: DO J1=1,NATOMS 
609:       IF (.NOT.ATOMACTIVE(J1)) CYCLE chicheck 
610:       IF (NBONDED(J1).NE.4) CYCLE chicheck 
611:       DO J2=1,4 
612:          IF (.NOT.ATOMACTIVE(BONDEDLIST(J1,J2))) CYCLE chicheck 
613:       ENDDO 
614: !     WRITE(*,'(A,I6,A)') 'intlbfgs> Active four-coordinate atom ',J1,' has four active neighbours' 
615:       CHIRALACTIVE(J1)=.TRUE. 
616:       DO J3=1,INTIMAGE+2 
617:          CENTRE_COORDS(1)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+1) 
618:          CENTRE_COORDS(2)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+2) 
619:          CENTRE_COORDS(3)=XYZ(3*NATOMS*(J3-1)+3*(J1-1)+3) 
620:  
621:          DO J4=1,4 
622:             J2=BONDEDLIST(J1,J4) 
623:             NEIGHBOUR_COORDS(3*(J4-1)+1)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+1) 
624:             NEIGHBOUR_COORDS(3*(J4-1)+2)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+2) 
625:             NEIGHBOUR_COORDS(3*(J4-1)+3)=XYZ(3*NATOMS*(J3-1)+3*(J2-1)+3) 
626:          ENDDO 
627:          CHIRALSR=CHIRALITY_SR(NEIGHBOUR_COORDS,CENTRE_COORDS) 
628: !        WRITE(*,'(A,I6,A,I6,A,L5)') 'intlbfgs> Atom ',J2,' image ',J3,' chirality=',CHIRALSR 
629:          IF ((J3.GT.1).AND.(CHIRALSR.NEQV.CHIRALSRP)) THEN 
630:             WRITE(*,'(A,I6,A,I6,A,I6)') 'intlbfgs> Atom ',J2,' image ',J3,' chirality CHANGED; use previous image coordinates'   
631:             DO J4=1,4 
632:                J2=BONDEDLIST(J1,J4) 
633:                XYZ(3*NATOMS*(J3-1)+3*(J2-1)+1)=XYZ(3*NATOMS*(J3-2)+3*(J2-1)+1) 
634:                XYZ(3*NATOMS*(J3-1)+3*(J2-1)+2)=XYZ(3*NATOMS*(J3-2)+3*(J2-1)+2) 
635:                XYZ(3*NATOMS*(J3-1)+3*(J2-1)+3)=XYZ(3*NATOMS*(J3-2)+3*(J2-1)+3) 
636:             ENDDO 
637:          ENDIF 
638:          IF (J3.EQ.1) CHIRALSRP=CHIRALSR  ! just use result for fixed end point image 1 
639:       ENDDO 
640:    ENDDO chicheck 
641:  
642: !  STOP 
643:  
644:    NDUMMY=1570:    NDUMMY=1
645:    DO J1=1,NPERMGROUP571:    DO J1=1,NPERMGROUP
646:       IF (GROUPACTIVE(J1)) GOTO 975572:       IF (GROUPACTIVE(J1)) GOTO 975
647:       DO J2=1,NPERMSIZE(J1)573:       DO J2=1,NPERMSIZE(J1)
648:          IF (.NOT.ATOMACTIVE(PERMGROUP(NDUMMY+J2-1))) GOTO 975574:          IF (.NOT.ATOMACTIVE(PERMGROUP(NDUMMY+J2-1))) GOTO 975
649:       ENDDO575:       ENDDO
650:       GROUPACTIVE(J1)=.TRUE.576:       GROUPACTIVE(J1)=.TRUE.
651:       IF (DEBUG) WRITE(*,'(A,I6,A)') ' intlbfgs> All permutable atoms in group ',J1,' are active'577:       WRITE(*,'(A,I6,A)') ' intlbfgs> All permutable atoms in group ',J1,' are active'
652: 975   NDUMMY=NDUMMY+NPERMSIZE(J1)578: 975   NDUMMY=NDUMMY+NPERMSIZE(J1)
653:    ENDDO 579:    ENDDO 
654: 580: 
655:    IF (NITERDONE.EQ.1) THEN581: !  COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint
656:       COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint582: !  COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint
657:       COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint583: !  WRITE(*,'(A)') ' intlbfgs> checking alignment of endpoints - should be no permutations'
658:       WRITE(*,'(A)') ' intlbfgs> checking alignment of endpoints - should be no permutations'584: !  CALL LOPERMDIST(COORDSB,COORDSC,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,.FALSE.,RMATBEST,0,NMOVE,PERM)
659:       CALL LOPERMDIST(COORDSB,COORDSC,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCE,DIST2,.FALSE.,RMATBEST,0,NMOVE,PERM)585: !  WRITE(*,'(A,G20.10,A,I6)') ' intlbfgs> endpoint distance ',DISTANCE,' permutations=',NMOVE 
660:       WRITE(*,'(A,G20.10,A,I6)') ' intlbfgs> endpoint distance ',DISTANCE,' permutations=',NMOVE  
661:    ENDIF 
662: 586: 
 587:    COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint
 588:    COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint
663:    SAVELOCALPERMCUT=LOCALPERMCUT589:    SAVELOCALPERMCUT=LOCALPERMCUT
 590: !
 591: ! needs to be a bit sloppier than usual to allow for distorted geometries. Otherwise we may get wrong assignments
 592: ! because we don't have enough neighbours.
 593: !
664:    LOCALPERMCUT=QCIPERMCUT ! 0.8D0 594:    LOCALPERMCUT=QCIPERMCUT ! 0.8D0 
665:    np: DO J1=1,NPERMGROUP595:    DO J3=1,INTIMAGE
666:       IF (.NOT.GROUPACTIVE(J1)) CYCLE np596:       np: DO J1=1,NPERMGROUP
667:       DO J3=1,INTIMAGE597:          IF (.NOT.GROUPACTIVE(J1)) CYCLE np
668: !        WRITE(*,'(A,I6,A,I6)') 'intlbfgs> Doing group ',J1,' image ',J3+1 
669:  
670: !          COORDSB(1:3*NATOMS)=XYZ(3*NATOMS*(J3-1)+1:3*NATOMS*J3) ! coordinates for intervening image J3-1, which is the starting endpoint for J3=1 
671: !          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))  ! coordinates for intervening image J3, which is image J3+1 including starting endpoint 
672: !          CALL QCIOPTPERM(COORDSB,COORDSA,NATOMS,DEBUG,J1,PERM,STARTGROUP,J3) 
673: !          IF (.NOT.IDENTITY) THEN 
674: !             DO J2=1,NATOMS 
675: ! !           IF (PERM(J2).NE.J2) WRITE(*,'(4(A,I6))') ' intlbfgs> image ',J3+1,' qcipermopt would move atom ',J2,' to position ',PERM(J2), &   
676: ! !  &   ' group ',J1   
677: !                IF (PERM(J2).NE.J2) WRITE(*,'(4(A,I6))') ' intlbfgs> image ',J3+1,' qcipermopt move atom ',J2,' to position ',PERM(J2), & 
678: !    &   ' group ',J1   
679: !                COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERM(J2)-1)+1) 
680: !                COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERM(J2)-1)+2) 
681: !                COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERM(J2)-1)+3) 
682: !             ENDDO 
683: !             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS) 
684: !          ELSE 
685: !             WRITE(*,'(A,I6,A,I6)') 'intlbfgs> identity permutation for image ',J3+1,' with preceeding image ',J3 
686: !          ENDIF 
687:  
688:          COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint 
689:          COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint 
690:          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint598:          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint
691:          CALL LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCES,DIST2,.FALSE.,RMATBEST,J1,NMOVES,PERMS)599:          CALL LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCES,DIST2,.FALSE.,RMATBEST,J1,NMOVES,PERMS)
692:          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint600:          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint
693:          CALL LOPERMDIST(COORDSC,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCEF,DIST2,.FALSE.,RMATBEST,J1,NMOVEF,PERMF)601:          CALL LOPERMDIST(COORDSC,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCEF,DIST2,.FALSE.,RMATBEST,J1,NMOVEF,PERMF)
694: !        WRITE(*,'(A,I6,A,I6,A,2G20.10,A,2I6)') ' intlbfgs> image ',J3+1,' group ',J1,' start and finish distances=', &602:          IF ((NMOVES.GT.0).AND.(NMOVEF.GT.0)) THEN
695: !  &                                       DISTANCES,DISTANCEF,' permutations=',NMOVES, NMOVEF603:             WRITE(*,'(A,I6,A,I6,A,2G20.10,A,2I6)') ' intlbfgs> image ',J3+1,' group ',J1,' start and finish distances =', &
696: 604:   &                                       DISTANCES,DISTANCEF,' permutations=',NMOVES, NMOVEF
697:          IF ((NMOVES.GT.0).OR.(NMOVEF.GT.0)) THEN 
698:             WRITE(*,'(A,I6,A,I6,A,2G20.10,A,2I6)') ' intlbfgs> image ',J3+1,' group ',J1,' start and finish distances=', & 
699:    &                                           DISTANCES,DISTANCEF,' permutations=',NMOVES, NMOVEF 
700:          ENDIF605:          ENDIF
701:  
702:          IF ((NMOVES.GT.0).AND.(NMOVES.EQ.NMOVEF)) THEN606:          IF ((NMOVES.GT.0).AND.(NMOVES.EQ.NMOVEF)) THEN
703:             COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) 607:             COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) 
704:             DO J2=1,NATOMS608:             DO J2=1,NATOMS
705:                IF (PERMS(J2).NE.J2) THEN609:                IF (PERMS(J2).NE.J2) THEN
706:                   IF (PERMF(J2).EQ.PERMS(J2)) THEN610:                   IF (PERMF(J2).EQ.PERMS(J2)) THEN
707:                      WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &611:                      WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &
708:    &                 ' consistent non-identity lopermdist permutation for start and finish, move atom ',PERMS(J2),' to position ',J2   612:   &                                           ' consistent non-identity permutation for start and finish, move atom ', &  
709: !                WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, &613:   &                                           PERMS(J2),' to position ',J2  
710: !  &           ' consistent non-identity lopermdist permutation for start and finish, would move atom ',PERMS(J2),' to position ',J2   
711:                      COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+1)614:                      COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+1)
712:                      COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+2)615:                      COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+2)
713:                      COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+3)616:                      COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+3)
714:                   ELSE617:                   ELSE
715:                      WRITE(*,'(A,I6,A,I6)') ' intlbfgs> inconsistent non-identity lopermdist permutations for start and finish'618:                      WRITE(*,'(A,I6,A,I6)') ' intlbfgs> inconsistent non-identity permutations for start and finish'
716:                   ENDIF619:                   ENDIF
717:                ENDIF620:                ENDIF
718:             ENDDO621:             ENDDO
719:             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS)622:             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS)
720:          ENDIF623:          ENDIF
721:        ENDDO624:       ENDDO np
722:     ENDDO np625:    ENDDO
723:     LOCALPERMCUT=SAVELOCALPERMCUT626:    LOCALPERMCUT=SAVELOCALPERMCUT
724: 627:    IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)
725: !STOP628: !  STOP
726: !  
727: !    COORDSB(1:3*NATOMS)=XYZ(1:3*NATOMS) ! starting endpoint 
728: !    COORDSC(1:3*NATOMS)=XYZ(3*NATOMS*(INTIMAGE+1)+1:3*NATOMS*(INTIMAGE+2)) ! finish endpoint 
729: !    SAVELOCALPERMCUT=LOCALPERMCUT 
730: ! ! 
731: ! ! needs to be a bit sloppier than usual to allow for distorted geometries. Otherwise we may get wrong assignments 
732: ! ! because we don't have enough neighbours. 
733: ! ! 
734: !    LOCALPERMCUT=QCIPERMCUT ! 0.8D0  
735: !    DO J3=1,INTIMAGE 
736: !       np: DO J1=1,NPERMGROUP 
737: !          IF (.NOT.GROUPACTIVE(J1)) CYCLE np 
738: !          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint 
739: !          CALL LOPERMDIST(COORDSB,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCES,DIST2,.FALSE.,RMATBEST,J1,NMOVES,PERMS) 
740: !          COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1)) ! coordinates for intervening image J3, which is image J3+1 including starting endpoint 
741: !          CALL LOPERMDIST(COORDSC,COORDSA,NATOMS,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT,TWOD,DISTANCEF,DIST2,.FALSE.,RMATBEST,J1,NMOVEF,PERMF) 
742: !             WRITE(*,'(A,I6,A,I6,A,2G20.10,A,2I6)') ' intlbfgs> image ',J3+1,' group ',J1,' start and finish distances=', & 
743: !   &                                       DISTANCES,DISTANCEF,' permutations=',NMOVES, NMOVEF 
744: !          IF ((NMOVES.GT.0).AND.(NMOVEF.GT.0)) THEN 
745: !             WRITE(*,'(A,I6,A,I6,A,2G20.10,A,2I6)') ' intlbfgs> image ',J3+1,' group ',J1,' start and finish distances=', & 
746: !   &                                       DISTANCES,DISTANCEF,' permutations=',NMOVES, NMOVEF 
747: !          ENDIF 
748: !          IF ((NMOVES.GT.0).AND.(NMOVES.EQ.NMOVEF)) THEN 
749: !             COORDSA(1:3*NATOMS)=XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))  
750: !             DO J2=1,NATOMS 
751: !                IF (PERMS(J2).NE.J2) THEN 
752: !                   IF (PERMF(J2).EQ.PERMS(J2)) THEN 
753: !                      WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> image ',J3+1, & 
754: !   &                                           ' consistent non-identity permutation for start and finish, move atom ', &   
755: !   &                                           PERMS(J2),' to position ',J2   
756: !                      COORDSA(3*(J2-1)+1)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+1) 
757: !                      COORDSA(3*(J2-1)+2)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+2) 
758: !                      COORDSA(3*(J2-1)+3)=XYZ(3*NATOMS*J3+3*(PERMS(J2)-1)+3) 
759: !                   ELSE 
760: !                      WRITE(*,'(A,I6,A,I6)') ' intlbfgs> inconsistent non-identity permutations for start and finish' 
761: !                   ENDIF 
762: !                ENDIF 
763: !             ENDDO 
764: !             XYZ(3*NATOMS*J3+1:3*NATOMS*(J3+1))=COORDSA(1:3*NATOMS) 
765: !          ENDIF 
766: !       ENDDO np 
767: !    ENDDO 
768: !    LOCALPERMCUT=SAVELOCALPERMCUT 
769: !    CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER) 
770: ! !  STOP 
771:  
772: ENDIF629: ENDIF
773: 630: 
 631: 
774: !632: !
775: !  Add next atom to active set if ADDATOM is true. 633: !  Add next atom to active set if ADDATOM is true. 
776: !  Constraints to atoms already in the active set are turned on634: !  Constraints to atoms already in the active set are turned on
777: !  and short-range repulsions to active atoms that are not distance constrained are turned on.635: !  and short-range repulsions to active atoms that are not distance constrained are turned on.
778: !  *** OLD Find nearest atom to active set attached by a constraint636: !  *** OLD Find nearest atom to active set attached by a constraint
779: !  *** NEW Find atom with most constraints to active set637: !  *** NEW Find atom with most constraints to active set
780: !  Turn on constraint terms for this atom with all previous members of the active set638: !  Turn on constraint terms for this atom with all previous members of the active set
781: !  Add repulsions to non-constrained atoms in this set639: !  Add repulsions to non-constrained atoms in this set
782: !  NTOADD is the number of atoms to add to the active set in each pass. 1 seems best!640: !  NTOADD is the number of atoms to add to the active set in each pass. 1 seems best!
783: !641: !
918:       IF (DEBUG) PRINT '(2(A,I6))', ' intlbfgs> Number of frozen images=',NIMAGEFREEZE,' / ',INTIMAGE776:       IF (DEBUG) PRINT '(2(A,I6))', ' intlbfgs> Number of frozen images=',NIMAGEFREEZE,' / ',INTIMAGE
919:    ENDIF777:    ENDIF
920:    !  We now have the proposed step - update geometry and calculate new gradient778:    !  We now have the proposed step - update geometry and calculate new gradient
921:    NDECREASE=0779:    NDECREASE=0
922: 20 X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)780: 20 X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
923: 781: 
924: !  IF (.NOT.SWITCHED) THEN782: !  IF (.NOT.SWITCHED) THEN
925:    IF (.TRUE.) THEN783:    IF (.TRUE.) THEN
926: !     IF ((RMS.LT.INTRMSTOL*1.0D10).AND.(MOD(NITERDONE,10).EQ.0).AND.(NSTEPSMAX-NITERDONE.GT.100)) &784: !     IF ((RMS.LT.INTRMSTOL*1.0D10).AND.(MOD(NITERDONE,10).EQ.0).AND.(NSTEPSMAX-NITERDONE.GT.100)) &
927: ! &               CALL CHECKSEP(NMAXINT,NMININT,INTIMAGE,XYZ,(3*NATOMS),NATOMS)785: ! &               CALL CHECKSEP(NMAXINT,NMININT,INTIMAGE,XYZ,(3*NATOMS),NATOMS)
928: !     PRINT '(A,3I10)','NITERDONE,INTIMAGECHECK,mod=',NITERDONE,INTIMAGECHECK,MOD(NITERDONE,INTIMAGECHECK) 
929:       IF (REMOVEIMAGE.OR.(MOD(NITERDONE,INTIMAGECHECK).EQ.0)) THEN786:       IF (REMOVEIMAGE.OR.(MOD(NITERDONE,INTIMAGECHECK).EQ.0)) THEN
930: 864      CONTINUE ! for adding more than one image at a time787: 864      CONTINUE ! for adding more than one image at a time
931:          DMAX=-1.0D0788:          DMAX=-1.0D0
932:          ADMAX=-1.0D0789:          ADMAX=-1.0D0
933:          DMIN=HUGE(1.0D0)790:          DMIN=HUGE(1.0D0)
934:          DO J1=1,INTIMAGE+1791:          DO J1=1,INTIMAGE+1
935:             DUMMY=0.0D0792:             DUMMY=0.0D0
936: !           DO J2=1,3*NATOMS793: !           DO J2=1,3*NATOMS
937: !              IF (ATOMACTIVE((J2-1)/3+1)) THEN794: !              IF (ATOMACTIVE((J2-1)/3+1)) THEN
938: !                 DUMMY=DUMMY+( XYZ((J1-1)*3*NATOMS+J2) - XYZ(J1*3*NATOMS+J2) )**2795: !                 DUMMY=DUMMY+( XYZ((J1-1)*3*NATOMS+J2) - XYZ(J1*3*NATOMS+J2) )**2
966: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1823: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1
967:          ENDDO824:          ENDDO
968: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &825: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &
969: ! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE826: ! &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE
970: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &827: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &
971: ! &                                                  DMAX,' for images ',JMAX,JMAX+1828: ! &                                                  DMAX,' for images ',JMAX,JMAX+1
972: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> smallest image separation is ', &829: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> smallest image separation is ', &
973: ! &                                                  DMIN,' for images ',JMIN,JMIN+1830: ! &                                                  DMIN,' for images ',JMIN,JMIN+1
974: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,G20.10)') ' intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)831: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,G20.10)') ' intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)
975: !        IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN832: !        IF ((DMAX.GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE)) THEN
976: !        PRINT '(A,2G20.10)','SQRT(ADMAX),IMSEPMAX=',SQRT(ADMAX),IMSEPMAX 
977:          IF ((.NOT.REMOVEIMAGE).AND.((SQRT(ADMAX).GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE))) THEN833:          IF ((.NOT.REMOVEIMAGE).AND.((SQRT(ADMAX).GT.IMSEPMAX).AND.(INTIMAGE.LT.MAXINTIMAGE))) THEN
978:             JMAX=JA1834:             JMAX=JA1
979:             WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> Add an image between ',JMAX,' and ',JMAX+1,' INTIMAGE=',INTIMAGE835:             WRITE(*,'(A,I6,A,I6,A,I6)') ' intlbfgs> Add an image between ',JMAX,' and ',JMAX+1,' INTIMAGE=',INTIMAGE
980:             NITERUSE=0836:             NITERUSE=0
981:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))837:             ALLOCATE(DPTMP(3*NATOMS*(INTIMAGE+2)))
982:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))838:             DPTMP(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
983:             DEALLOCATE(XYZ)839:             DEALLOCATE(XYZ)
984:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))840:             ALLOCATE(XYZ(3*NATOMS*(INTIMAGE+3)))
985:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)841:             XYZ(1:3*NATOMS*JMAX)=DPTMP(1:3*NATOMS*JMAX)
986:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &842:             XYZ(3*NATOMS*JMAX+1:3*NATOMS*(JMAX+1))=(DPTMP(3*NATOMS*(JMAX-1)+1:3*NATOMS*JMAX) &
1175:                DMIN=DUMMY1031:                DMIN=DUMMY
1176:                JMIN=J11032:                JMIN=J1
1177:             ENDIF1033:             ENDIF
1178: !            IF (DEBUG) WRITE(*,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &1034: !            IF (DEBUG) WRITE(*,'(A,I6,A,I6,A,G20.10)')' intlbfgs> distance between images ', &
1179: !  &                                                  J1,' and ',J1+1,' is ',DUMMY1035: !  &                                                  J1,' and ',J1+1,' is ',DUMMY
1180: !            IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &1036: !            IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6)')' intlbfgs> largest atomic distance between images so far is ', &
1181: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+11037: !  &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1
1182:          ENDDO1038:          ENDDO
1183:          IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &1039:          IF (DEBUG) WRITE(*,'(A,G20.10,A,I6,A,2I6,A,I6)')' intlbfgs> largest atomic distance between images is ', &
1184:   &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE1040:   &                                                  SQRT(ADMAX),' for atom ',JA2,' and images ',JA1,JA1+1,' total images=',INTIMAGE
1185: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &1041:          IF (DEBUG) WRITE(*,'(A,G20.10,A,2I6)')' intlbfgs> largest image separation is ', &
1186: ! &                                                  DMAX,' for images ',JMAX,JMAX+11042:   &                                                  DMAX,' for images ',JMAX,JMAX+1
1187: !        IF (DEBUG) WRITE(*,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)1043:          IF (DEBUG) WRITE(*,'(A,G20.10,A,G20.10)') 'intlbfgs> Mean image separation=',DUMMY2/(INTIMAGE+1),' per active atom=',DUMMY2/((INTIMAGE+1)*NACTIVE)
1188: !        IF (SQRT(ADMAX).GT.IMSEPMAX) THEN1044: !        IF (SQRT(ADMAX).GT.IMSEPMAX) THEN
1189: !           KINT=MIN(1.0D6,KINT*1.1D0)1045: !           KINT=MIN(1.0D6,KINT*1.1D0)
1190: !        ELSE1046: !        ELSE
1191: !           KINT=MAX(1.0D-6,KINT/1.1D0)1047: !           KINT=MAX(1.0D-6,KINT/1.1D0)
1192: !        ENDIF1048: !        ENDIF
1193: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT1049: !        WRITE(*,'(A,G20.10)') 'intlbfgs> Spring constant is now ',KINT
1194:       ENDIF1050:       ENDIF
1195:    ENDIF1051:    ENDIF
1196: !1052: !
1197: ! End of add/subtract images block.1053: ! End of add/subtract images block.
1198: !1054: !
1199: !1055: !
1200: ! The new QCILPERMDIST check should be used instead of these lines1056: ! The new QCILPERMDIST check should be used instead of these lines
1201: !1057: !
1202: !  IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN1058: !  IF (QCIPERMCHECK.AND.(MOD(NITERDONE,QCIPERMCHECKINT).EQ.0)) THEN
1203: !     LDEBUG=.FALSE.1059: !     LDEBUG=.FALSE.
1204: !     DO J2=2,INTIMAGE+21060: !     DO J2=2,INTIMAGE+2
1205: !        CALL MINPERMDIST(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &1061: !        CALL ALIGN_DECIDE(XYZ((J2-2)*3*NATOMS+1:(J2-1)*3*NATOMS),XYZ((J2-1)*3*NATOMS+1:J2*3*NATOMS),NATOMS,LDEBUG, &
1206: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)1062: ! &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,DIST,DIST2,RIGIDBODY,RMAT)
1207: !     ENDDO1063: !     ENDDO
1208: !  ENDIF1064: !  ENDIF
1209: 1065: 
1210:    IF (.NOT.SWITCHED) THEN1066:    IF (.NOT.SWITCHED) THEN
1211:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)1067:       IF (MOD(NITERDONE,CHECKREPINTERVAL).EQ.0) CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),0,1)
1212:       IF (QCIADDREP.GT.0) THEN1068:       IF (QCIADDREP.GT.0) THEN
1213:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1069:          CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1214:       ELSEIF (CHECKCONINT) THEN1070:       ELSEIF (CHECKCONINT) THEN
1215:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)1071:          CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
1309:          MAXRMS=DUMMY1165:          MAXRMS=DUMMY
1310:          JMAXRMS=J11166:          JMAXRMS=J1
1311:       ENDIF1167:       ENDIF
1312:    ENDDO1168:    ENDDO
1313:    MAXRMS=SQRT(MAXRMS/(3*NACTIVE))1169:    MAXRMS=SQRT(MAXRMS/(3*NACTIVE))
1314:    SUMEEE=SUMEEE/INTIMAGE1170:    SUMEEE=SUMEEE/INTIMAGE
1315:    SUMEEE2=SUMEEE2+EEE(J1)**21171:    SUMEEE2=SUMEEE2+EEE(J1)**2
1316:    SUMEEE2=SUMEEE2/INTIMAGE1172:    SUMEEE2=SUMEEE2/INTIMAGE
1317:    SIGMAEEE=SQRT(SUMEEE2-SUMEEE**2)1173:    SIGMAEEE=SQRT(SUMEEE2-SUMEEE**2)
1318:    REMOVEIMAGE=.FALSE.1174:    REMOVEIMAGE=.FALSE.
1319: !    IF (ABS(MAXEEE-SUMEEE).GT.3.0D0*SIGMAEEE) THEN1175: !  IF (ABS(MAXEEE-SUMEEE).GT.3.0D0*SIGMAEEE) THEN
1320: ! !  IF (MAXEEE.GT.1.0D2*MINEEE) THEN1176: !  IF (MAXEEE.GT.1.0D2*MINEEE) THEN
1321:       WRITE(*,'(A,I8,A,G20.10,A,G20.10,A)') ' intlbfgs> Highest image ',JMAXEEE,' energy ',MAXEEE,' is ',ABS(MAXEEE-SUMEEE)/SIGMAEEE, &1177:       WRITE(*,'(A,I8,A,G20.10,A,G20.10,A)') ' intlbfgs> Highest image ',JMAXEEE,' energy ',MAXEEE,' is ',ABS(MAXEEE-SUMEEE)/SIGMAEEE, &
1322:   &                        ' sigma from the mean'1178:   &                        ' sigma from the mean'
1323: !       REMOVEIMAGE=.TRUE.1179: !     REMOVEIMAGE=.TRUE.
1324: !       IF (NITERDONE-NLASTGOODE.LT.20) REMOVEIMAGE=.FALSE.1180: !     IF (NITERDONE-NLASTGOODE.LT.40) REMOVEIMAGE=.FALSE.
1325: !    ENDIF1181: !  ENDIF
1326: 1182: 
1327:    IF (DEBUG) THEN1183:    IF (DEBUG) THEN
1328: !     WRITE(*,'(A,2G20.10)') ' intlbfgs> mean and sigma of image energies=',SUMEEE,SIGMAEEE1184:       WRITE(*,'(A,2G20.10)') ' intlbfgs> mean and sigma of image energies=',SUMEEE,SIGMAEEE
1329:       WRITE(*,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &1185:       WRITE(*,'(A,I6,2G20.10,3(G20.10,I8))') ' intlbfgs> steps: ',NITERDONE,ETOTAL/INTIMAGE,RMS,STEPTOT,NACTIVE, &
1330:   &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS1186:   &                                                        MAXEEE,JMAXEEE,MAXRMS,JMAXRMS
1331:       CALL FLUSH(6)1187:       CALL FLUSH(6)
1332:    ENDIF1188:    ENDIF
1333: 1189: 
1334:    IF (.NOT.SWITCHED) THEN1190:    IF (.NOT.SWITCHED) THEN
1335: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN1191: !     IF ((NITERDONE-NLASTGOODE.GT.INTRELSTEPS).AND.((ETOTAL.GT.LASTGOODE).OR.(ETOTAL/INTIMAGE.GT.MAXCONE*1.0D8))) THEN
1336:       IF (.FALSE.) THEN ! no backtracking1192:       IF (.FALSE.) THEN ! no backtracking
1337:          WRITE(*,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE1193:          WRITE(*,'(2(A,I6))') ' intlbfgs> Backtracking ',NBACKTRACK,' steps, current active atoms=',NACTIVE
1338:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+11194:          NTRIES(NEWATOM)=NTRIES(NEWATOM)+1
1443:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1299:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1444:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 1300:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 
1445:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 1301:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 
1446:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21302:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1447:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1303:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
1448: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX1304: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX
1449: 1305: 
1450:    IF (EXITSTATUS > 0) THEN  1306:    IF (EXITSTATUS > 0) THEN  
1451:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1307:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1452: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771308: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1453:          PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.2D0,NACTIVE*1.0D0/(NATOMS*1.0D0))1309:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1454:          IF (MAXEEE.GT.MAXCONE*MAX(0.2D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777 
1455:          IF (NACTIVE.LT.NATOMS) THEN 1310:          IF (NACTIVE.LT.NATOMS) THEN 
1456:             ADDATOM=.TRUE.1311:             ADDATOM=.TRUE.
1457:             GOTO 7771312:             GOTO 777
1458:          ENDIF1313:          ENDIF
1459:          CALL MYCPU_TIME(FTIME,.FALSE.)1314:          CALL MYCPU_TIME(FTIME,.FALSE.)
1460:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1315:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,F10.1)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1461:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1316:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1462:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)1317:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER)
1463:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1318:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1464:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1319:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
1465:          DO J1=1,NATOMS1320:          DO J1=1,NATOMS
1466:             IF (.NOT.ATOMACTIVE(J1)) THEN1321:             IF (.NOT.ATOMACTIVE(J1)) THEN
1467:                WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> ERROR *** number of active atoms=',NACTIVE,' but atom ',J1,' is not active'1322:                WRITE(*,'(A,I6,A,I6,A)') ' intlbfgs> ERROR *** number of active atoms=',NACTIVE,' but atom ',J1,' is not active'
1468:             ENDIF1323:             ENDIF
1469:          ENDDO1324:          ENDDO
1470:          NSTEPSMAX=NITERDONE+INTCONSTEPS1325:          NSTEPSMAX=NITERDONE+INTCONSTEPS
1581: 1436: 
1582: NNREPULSIVE=NNSTART1437: NNREPULSIVE=NNSTART
1583: DO JJ=NSTART,NREPULSIVE1438: DO JJ=NSTART,NREPULSIVE
1584:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**21439:    COMPARE=(CHECKREPCUTOFF*REPCUT(JJ))**2
1585:    NI=REPI(JJ)1440:    NI=REPI(JJ)
1586:    NJ=REPJ(JJ)1441:    NJ=REPJ(JJ)
1587:    DO KK=1,INTIMAGE+2 ! first check for standard distances within threshold1442:    DO KK=1,INTIMAGE+2 ! first check for standard distances within threshold
1588:       LDIST=(XYZ((KK-1)*NOPT+3*(NI-1)+1)-XYZ((KK-1)*NOPT+3*(NJ-1)+1))**2 &1443:       LDIST=(XYZ((KK-1)*NOPT+3*(NI-1)+1)-XYZ((KK-1)*NOPT+3*(NJ-1)+1))**2 &
1589:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+2)-XYZ((KK-1)*NOPT+3*(NJ-1)+2))**2 &1444:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+2)-XYZ((KK-1)*NOPT+3*(NJ-1)+2))**2 &
1590:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+3)-XYZ((KK-1)*NOPT+3*(NJ-1)+3))**21445:   &        +(XYZ((KK-1)*NOPT+3*(NI-1)+3)-XYZ((KK-1)*NOPT+3*(NJ-1)+3))**2
1591: !     IF ((REPI(JJ).EQ.2024).OR.(REPJ(JJ).EQ.2024)) THEN 
1592: !        WRITE(*,'(A,2I12,2I6,3G20.10)') 'checkrep> KK,JJ,repi,repj,ldist,compare,repcut=', & 
1593: ! &                                    KK,JJ,repi(JJ),repj(JJ),ldist,compare,repcut(jj) 
1594: !     ENDIF 
1595:       IF (LDIST.LT.COMPARE) THEN1446:       IF (LDIST.LT.COMPARE) THEN
1596:          NNREPULSIVE=NNREPULSIVE+11447:          NNREPULSIVE=NNREPULSIVE+1
1597:          NREPI(NNREPULSIVE)=NI1448:          NREPI(NNREPULSIVE)=NI
1598:          NREPJ(NNREPULSIVE)=NJ1449:          NREPJ(NNREPULSIVE)=NJ
1599:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)1450:          NREPCUT(NNREPULSIVE)=REPCUT(JJ)
1600: !        IF ((REPI(JJ).EQ.2024).OR.(REPJ(JJ).EQ.2024)) THEN 
1601: !           WRITE(*,'(A)') 'checkrep> KK,JJ,ldist < compare : active' 
1602: !        ENDIF 
1603:          GOTO 2461451:          GOTO 246
1604:       ENDIF1452:       ENDIF
1605:    ENDDO 1453:    ENDDO 
1606:    COMPARE=CHECKREPCUTOFF*REPCUT(JJ)1454:    COMPARE=CHECKREPCUTOFF*REPCUT(JJ)
1607:    DO KK=2,INTIMAGE+2 ! now check internal minima within threshold1455:    DO KK=2,INTIMAGE+2 ! now check internal minima within threshold
1608:       DMIN=1.0D101456:       DMIN=1.0D10
1609:       NI2=NOPT*(KK-2)+3*(NI-1)1457:       NI2=NOPT*(KK-2)+3*(NI-1)
1610:       NI1=NOPT*(KK-1)+3*(NI-1)1458:       NI1=NOPT*(KK-1)+3*(NI-1)
1611:       NJ2=NOPT*(KK-2)+3*(NJ-1)1459:       NJ2=NOPT*(KK-2)+3*(NJ-1)
1612:       NJ1=NOPT*(KK-1)+3*(NJ-1)1460:       NJ1=NOPT*(KK-1)+3*(NJ-1)
1731: WRITE(UNIT=LUNIT,FMT='(2G24.13)') EEE(INTIMAGE+2)1579: WRITE(UNIT=LUNIT,FMT='(2G24.13)') EEE(INTIMAGE+2)
1732: 1580: 
1733: CLOSE(LUNIT)1581: CLOSE(LUNIT)
1734: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'1582: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'
1735: 1583: 
1736: END SUBROUTINE WRITEPROFILE1584: END SUBROUTINE WRITEPROFILE
1737: 1585: 
1738: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)1586: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE)
1739: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &1587: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &
1740:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &1588:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &
1741:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, ATOMSTORES, QCIADDACIDT, &1589:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, &
1742:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &1590:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &
1743:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, QCIADDREP1591:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, QCIADDREP
1744: USE COMMONS, ONLY: NATOMS, DEBUG1592: USE COMMONS, ONLY: NATOMS, DEBUG
1745: IMPLICIT NONE1593: IMPLICIT NONE
1746: INTEGER INTIMAGE1594: INTEGER INTIMAGE
1747: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &1595: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &
1748:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE1596:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE
1749: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN1597: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN
1750: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID1598: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS)
1751: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)1599: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)
1752: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID1600: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS)
1753: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE01601: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE0
1754: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE1602: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE
1755: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), XSAVED(3,INTIMAGE+2), XSAVEC(3,INTIMAGE+2), XSAVE0(3,INTIMAGE+2),FRAC,RAN1, &1603: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), XSAVED(3,INTIMAGE+2), XSAVEC(3,INTIMAGE+2), XSAVE0(3,INTIMAGE+2),FRAC,RAN1, &
1756:   &              RMS,EEE(INTIMAGE+2),GGG((3*NATOMS)*(INTIMAGE+2)),ETOTAL,DS,DF,DNORM1604:   &              RMS,EEE(INTIMAGE+2),GGG((3*NATOMS)*(INTIMAGE+2)),ETOTAL,DS,DF,DNORM
1757: 1605: 
1758:  
1759: NTOADD=11606: NTOADD=1
1760: NADDED=01607: NADDED=0
1761: CHOSENACID=.FALSE. 
1762:  
1763: ! 
1764: ! ATOMSTORES(J1) is the residue for atom J1 
1765: ! 
1766: ! AMBER12_GET_RESDATA needs a data type in amber12 interface and the number of residues 
1767: ! defined in amber12interface.f90 
1768: ! 
1769: 1608: 
1770: !1609: !
1771: ! Save current number of repulsions and number that are active to speed up the1610: ! Save current number of repulsions and number that are active to speed up the
1772: ! calls to CHECKREP1611: ! calls to CHECKREP
1773: !1612: !
1774: NNREPSAVE=NNREPULSIVE1613: NNREPSAVE=NNREPULSIVE
1775: NREPSAVE=NREPULSIVE1614: NREPSAVE=NREPULSIVE
1776: 542   CONTINUE1615: 542   CONTINUE
1777: !     DUMMY=1.0D1001616: !     DUMMY=1.0D100
1778:       NBEST=01617:       NBEST=0
1841: !          ENDIF1680: !          ENDIF
1842: !       ENDDO choosenew1681: !       ENDDO choosenew
1843: 1682: 
1844: !1683: !
1845: !  Choose NEWATOM deterministically. Take the inactive atom with the shortest constrained distance.1684: !  Choose NEWATOM deterministically. Take the inactive atom with the shortest constrained distance.
1846: !1685: !
1847:       DUMMY2=1.0D1001686:       DUMMY2=1.0D100
1848:       DO J1=1,NCONSTRAINT1687:       DO J1=1,NCONSTRAINT
1849:          IF (CONACTIVE(J1)) CYCLE1688:          IF (CONACTIVE(J1)) CYCLE
1850:          IF (ATOMACTIVE(CONJ(J1))) THEN1689:          IF (ATOMACTIVE(CONJ(J1))) THEN
1851:             IF (CHOSENACID.AND.(.NOT.(ATOMSTORES(CONI(J1)).EQ.ACID))) THEN1690:             IF (.NOT.ATOMACTIVE(CONI(J1))) THEN
1852:             ELSE1691:                IF (CONDISTREF(J1).LT.DUMMY2) THEN
1853:                IF (.NOT.ATOMACTIVE(CONI(J1))) THEN1692:                   DUMMY2=CONDISTREF(J1)
1854:                   IF (CONDISTREF(J1).LT.DUMMY2) THEN1693:                   NEWATOM=CONI(J1)
1855:                      DUMMY2=CONDISTREF(J1) 
1856:                      NEWATOM=CONI(J1) 
1857:                   ENDIF 
1858:                ENDIF1694:                ENDIF
1859:             ENDIF1695:             ENDIF
1860:          ELSEIF (ATOMACTIVE(CONI(J1))) THEN1696:          ELSEIF (ATOMACTIVE(CONI(J1))) THEN
1861:             IF (CHOSENACID.AND.(.NOT.(ATOMSTORES(CONJ(J1)).EQ.ACID))) THEN1697:             IF (.NOT.ATOMACTIVE(CONJ(J1))) THEN
1862:             ELSE1698:                IF (CONDISTREF(J1).LT.DUMMY2) THEN
1863:                IF (.NOT.ATOMACTIVE(CONJ(J1))) THEN1699:                   DUMMY2=CONDISTREF(J1)
1864:                   IF (CONDISTREF(J1).LT.DUMMY2) THEN1700:                   NEWATOM=CONJ(J1)
1865:                      DUMMY2=CONDISTREF(J1) 
1866:                      NEWATOM=CONJ(J1) 
1867:                   ENDIF 
1868:                ENDIF1701:                ENDIF
1869:             ENDIF1702:             ENDIF
1870:          ENDIF1703:          ENDIF
1871:       ENDDO1704:       ENDDO
1872:       IF (DEBUG) WRITE(*,'(3(A,I6),A,F15.5)') ' intlbfgs> Choosing new active atom ',NEWATOM,' new constraints=', &1705:       IF (DEBUG) WRITE(*,'(3(A,I6),A,F15.5)') ' intlbfgs> Choosing new active atom ',NEWATOM,' new constraints=', &
1873:   &                                       NCONTOACTIVE(NEWATOM),' maximum=',NBEST,' shortest constraint=',DUMMY21706:   &                                       NCONTOACTIVE(NEWATOM),' maximum=',NBEST,' shortest constraint=',DUMMY2
1874: !     IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,ETOTAL)1707: !     IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,ETOTAL)
1875: !     IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1708: !     IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1876:       IF (QCIADDACIDT.AND.(.NOT.CHOSENACID)) THEN1709: 
1877:          ACID=ATOMSTORES(NEWATOM) 
1878:          CHOSENACID=.TRUE. 
1879:       ENDIF 
1880:           1710:           
1881:       IF (NEWATOM*NBEST.EQ.0) THEN ! sanity check1711:       IF (NEWATOM*NBEST.EQ.0) THEN ! sanity check
1882:          WRITE(*,'(A,I6,A,2I6)') ' intlbfgs> ERROR *** new active atom not set'1712:          WRITE(*,'(A,I6,A,2I6)') ' intlbfgs> ERROR *** new active atom not set'
1883:          STOP1713:          STOP
1884:       ELSE1714:       ELSE
1885: !1715: !
1886: !  We need a sorted list of up to 3 active atoms, sorted according to how well the1716: !  We need a sorted list of up to 3 active atoms, sorted according to how well the
1887: !  end point distance is preserved, even if they don't satisfy the constraint 1717: !  end point distance is preserved, even if they don't satisfy the constraint 
1888: !  condition. We want three atoms to use for a local axis system in the interpolation.1718: !  condition. We want three atoms to use for a local axis system in the interpolation.
1889: !1719: !
2443:          ELSE 2273:          ELSE 
2444:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest constraints'2274:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest constraints'
2445:             DO J1=2,INTIMAGE+12275:             DO J1=2,INTIMAGE+1
2446:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVE0(1:3,J1)2276:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVE0(1:3,J1)
2447:             ENDDO2277:             ENDDO
2448:             ETOTAL=ESAVE02278:             ETOTAL=ESAVE0
2449:          ENDIF2279:          ENDIF
2450:       ENDIF2280:       ENDIF
2451:       NADDED=NADDED+12281:       NADDED=NADDED+1
2452:       IF (NADDED.LT.NTOADD) GOTO 5422282:       IF (NADDED.LT.NTOADD) GOTO 542
2453: ! 
2454: ! Check whether we've added all atoms in the amino acid corresponding to the new atom. If not, go back to the top 
2455: ! and choose the next candidate. 
2456: ! 
2457:       IF (QCIADDACIDT) THEN 
2458:          DO J1=1,NATOMS 
2459:             IF ((ATOMSTORES(J1).EQ.ACID).AND.(.NOT.(ATOMACTIVE(J1)))) GOTO 542 
2460:          ENDDO 
2461:          WRITE(*,'(A,I6,A)') 'doaddatom> All atoms of residue ',ACID,' are active' 
2462:       ENDIF 
2463: 2283: 
2464:       IF (QCIRADSHIFTT) THEN2284:       IF (QCIRADSHIFTT) THEN
2465:          WRITE(*,'(A,F15.5)') ' intlbfgs> Applying radial shift for unconstrained atoms of ',QCIRADSHIFT2285:          WRITE(*,'(A,F15.5)') ' intlbfgs> Applying radial shift for unconstrained atoms of ',QCIRADSHIFT
2466:          WRITE(*,'(20I6)') CONLIST(1:NCONFORNEWATOM)2286:          WRITE(*,'(20I6)') CONLIST(1:NCONFORNEWATOM)
2467:          DO J1=2,INTIMAGE+12287:          DO J1=2,INTIMAGE+1
2468:             scaleloop: DO J2=1,NATOMS2288:             scaleloop: DO J2=1,NATOMS
2469:                IF (.NOT.ATOMACTIVE(J2)) CYCLE scaleloop2289:                IF (.NOT.ATOMACTIVE(J2)) CYCLE scaleloop
2470:                IF (J2.EQ.NEWATOM) CYCLE scaleloop2290:                IF (J2.EQ.NEWATOM) CYCLE scaleloop
2471:                DO J3=1,NCONFORNEWATOM2291:                DO J3=1,NCONFORNEWATOM
2472:                   IF (CONLIST(J3).EQ.J2) CYCLE scaleloop2292:                   IF (CONLIST(J3).EQ.J2) CYCLE scaleloop
2506: END SUBROUTINE DOADDATOM2326: END SUBROUTINE DOADDATOM
2507: 2327: 
2508: SUBROUTINE CHECKPERC(LXYZ,LINTCONSTRAINTTOL,NQCIFREEZE,NCPFIT)2328: SUBROUTINE CHECKPERC(LXYZ,LINTCONSTRAINTTOL,NQCIFREEZE,NCPFIT)
2509: USE KEY, ONLY : ATOMACTIVE, NCONSTRAINT, INTFROZEN, CONI, CONJ, CONDISTREF, INTCONMAX, INTCONSTRAINTTOL, &2329: USE KEY, ONLY : ATOMACTIVE, NCONSTRAINT, INTFROZEN, CONI, CONJ, CONDISTREF, INTCONMAX, INTCONSTRAINTTOL, &
2510:   &             INTCONSEP, NCONGEOM, CONGEOM, CONIFIX, CONJFIX, CONDISTREFFIX, INTCONCUT, &2330:   &             INTCONSEP, NCONGEOM, CONGEOM, CONIFIX, CONJFIX, CONDISTREFFIX, INTCONCUT, &
2511:   &             NCONSTRAINTFIX, BULKT, TWOD, RIGIDBODY, CONDATT, CONCUT, CONCUTFIX, &2331:   &             NCONSTRAINTFIX, BULKT, TWOD, RIGIDBODY, CONDATT, CONCUT, CONCUTFIX, &
2512:   &             BONDS, QCIAMBERT, QCIADDREP, QCIADDREPCUT, QCIBONDS, QCISECOND2332:   &             BONDS, QCIAMBERT, QCIADDREP, QCIADDREPCUT, QCIBONDS, QCISECOND
2513: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM32333: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3
2514: IMPLICIT NONE2334: IMPLICIT NONE
2515: INTEGER NDIST1(NATOMS), NCYCLE, DMIN1, DMAX1, NUNCON1, J1, J2, J3, NQCIFREEZE, J4, NCPFIT, LUNIT, GETUNIT2335: INTEGER NDIST1(NATOMS), NCYCLE, DMIN1, DMAX1, NUNCON1, J1, J2, J3, NQCIFREEZE, J4, NCPFIT, LUNIT, GETUNIT
2516: INTEGER NI1, NJ1, NI2, NJ2, J5, ATOM1, ATOM2, ACID2336: INTEGER NI1, NJ1, NI2, NJ2, J5, ATOM1, ATOM2
2517: DOUBLE PRECISION LINTCONSTRAINTTOL, MAXCONDIST, MINCONDIST, DS, DF, LXYZ((3*NATOMS)*2)2337: DOUBLE PRECISION LINTCONSTRAINTTOL, MAXCONDIST, MINCONDIST, DS, DF, LXYZ((3*NATOMS)*2)
2518: DOUBLE PRECISION DSMIN, DSMAX, DSMEAN, D, DIST2, RMAT(3,3), DUMMY, X1, Y1, Z1, X2, Y2, Z2, DMIN, D22338: DOUBLE PRECISION DSMIN, DSMAX, DSMEAN, D, DIST2, RMAT(3,3), DUMMY, X1, Y1, Z1, X2, Y2, Z2, DMIN, D2
2519: LOGICAL CHANGED, LDEBUG, CONFILET2339: LOGICAL CHANGED, LDEBUG, CONFILET
2520: LOGICAL :: CALLED=.FALSE.2340: LOGICAL :: CALLED=.FALSE.
2521: SAVE CALLED2341: SAVE CALLED
2522: !for QCIAMBER2342: !for QCIAMBER
2523: INTEGER NBOND, NDUMMY2343: INTEGER NBOND, NDUMMY
2524: 2344: 
2525: LINTCONSTRAINTTOL=INTCONSTRAINTTOL2345: LINTCONSTRAINTTOL=INTCONSTRAINTTOL
2526: 2346: 
2565:       ENDDO2385:       ENDDO
2566:       NCONSTRAINT=J22386:       NCONSTRAINT=J2
2567:       WRITE(*,'(A,I6,A)') ' checkperc> After allowing for frozen atoms there are ',NCONSTRAINT,' constraints'2387:       WRITE(*,'(A,I6,A)') ' checkperc> After allowing for frozen atoms there are ',NCONSTRAINT,' constraints'
2568:       RETURN 2388:       RETURN 
2569:    ELSE2389:    ELSE
2570: !2390: !
2571: ! Put reference minima in optimal permutational alignment with reference minimum one.2391: ! Put reference minima in optimal permutational alignment with reference minimum one.
2572: !2392: !
2573:       DO J2=2,NCONGEOM2393:       DO J2=2,NCONGEOM
2574:          LDEBUG=.FALSE.2394:          LDEBUG=.FALSE.
2575:          CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),CONGEOM(J2,1:3*NATOMS),NATOMS,LDEBUG, &2395:          CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),CONGEOM(J2,1:3*NATOMS),NATOMS,LDEBUG, &
2576:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)2396:   &                       PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,DIST2,RIGIDBODY,RMAT)
2577:       ENDDO2397:       ENDDO
2578:    ENDIF2398:    ENDIF
2579:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))2399:    ALLOCATE(CONIFIX(INTCONMAX),CONJFIX(INTCONMAX),CONCUTFIX(INTCONMAX),CONDISTREFFIX(INTCONMAX))
2580: ENDIF2400: ENDIF
2581: 2401: 
2582: INQUIRE(FILE='constraintfile',EXIST=CONFILET)2402: INQUIRE(FILE='constraintfile',EXIST=CONFILET)
2583: 2403: 
2584: 51   NCONSTRAINT=0 2404: 51   NCONSTRAINT=0 
2585: MAXCONDIST=-1.0D02405: MAXCONDIST=-1.0D0
2667: !2487: !
2668:       DO2488:       DO
2669:          READ(LUNIT,*,END=534)  J2, J32489:          READ(LUNIT,*,END=534)  J2, J3
2670: !2490: !
2671: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to2491: ! Forbid constraints corresponding to atoms distant in sequence. Set INTCONSEP to number of sites to
2672: ! turn this off2492: ! turn this off
2673: !2493: !
2674:          IF (J3-J2.GT.INTCONSEP) CYCLE2494:          IF (J3-J2.GT.INTCONSEP) CYCLE
2675:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms2495:          IF (INTFROZEN(J2).AND.INTFROZEN(J3)) CYCLE ! no constraints between intfrozen atoms
2676:          NCONSTRAINT=NCONSTRAINT+12496:          NCONSTRAINT=NCONSTRAINT+1
 2497:          IF (DEBUG) WRITE(*,'(A,2I6,A,I6)') 'intlbfgs> Adding extra constraint for atoms ',J2,J3,'  total=',NCONSTRAINT
2677:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &2498:          DS=SQRT((LXYZ(3*(J2-1)+1)-LXYZ(3*(J3-1)+1))**2 &
2678:   &             +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 &2499:   &             +(LXYZ(3*(J2-1)+2)-LXYZ(3*(J3-1)+2))**2 &
2679:   &             +(LXYZ(3*(J2-1)+3)-LXYZ(3*(J3-1)+3))**2)2500:   &             +(LXYZ(3*(J2-1)+3)-LXYZ(3*(J3-1)+3))**2)
2680:          DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 &2501:          DF=SQRT((LXYZ(3*NATOMS+3*(J2-1)+1)-LXYZ(3*NATOMS+3*(J3-1)+1))**2 &
2681:   &             +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 &2502:   &             +(LXYZ(3*NATOMS+3*(J2-1)+2)-LXYZ(3*NATOMS+3*(J3-1)+2))**2 &
2682:   &             +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2)2503:   &             +(LXYZ(3*NATOMS+3*(J2-1)+3)-LXYZ(3*NATOMS+3*(J3-1)+3))**2)
2683:          IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE2504:          IF (NCONSTRAINT.GT.INTCONMAX) CALL CONDOUBLE
2684:          CONI(NCONSTRAINT)=J22505:          CONI(NCONSTRAINT)=J2
2685:          CONJ(NCONSTRAINT)=J32506:          CONJ(NCONSTRAINT)=J3
2686:          CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D02507:          CONDISTREF(NCONSTRAINT)=(DF+DS)/2.0D0
2687:          CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D02508:          CONCUT(NCONSTRAINT)=ABS(DF-DS)/2.0D0
2688:          IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)2509:          IF (CONDISTREF(NCONSTRAINT).GT.MAXCONDIST) MAXCONDIST=CONDISTREF(NCONSTRAINT)
2689:          IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)2510:          IF (CONDISTREF(NCONSTRAINT).LT.MINCONDIST) MINCONDIST=CONDISTREF(NCONSTRAINT)
2690:          IF (DEBUG) WRITE(*,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> extra constraint distance for ',CONI(NCONSTRAINT), &2511:          IF (DEBUG) WRITE(*,'(A,2I6,A,2F12.2,A,F12.4,A,I8)') ' intlbfgs> constrain distance for ',CONI(NCONSTRAINT), &
2691:   &                     CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &2512:   &                     CONJ(NCONSTRAINT),' values are ',DS,DF,' fraction=',2*ABS(DS-DF)/(DS+DF), &
2692:   &                  ' # constraints=',NCONSTRAINT2513:   &                  ' # constraints=',NCONSTRAINT
2693:       ENDDO2514:       ENDDO
2694: 534   CONTINUE2515: 534   CONTINUE
2695:       CLOSE(LUNIT)2516:       CLOSE(LUNIT)
2696:       IF (NCONSTRAINT-NDUMMY.GT.0) WRITE(*,'(A,I6,2(A,F15.5))') ' intlbfgs> Extra distance constraints: ',NCONSTRAINT-NDUMMY2517:       IF (NCONSTRAINT-NDUMMY.GT.0) WRITE(*,'(A,I6,2(A,F15.5))') ' intlbfgs> Extra distance constraints: ',NCONSTRAINT-NDUMMY
2697:       WRITE(*,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT,' shortest=',MINCONDIST,' longest=',MAXCONDIST2518:       WRITE(*,'(A,I6,2(A,F15.5))') ' intlbfgs> Total distance constraints=',NCONSTRAINT,' shortest=',MINCONDIST,' longest=',MAXCONDIST
2698:       CLOSE(LUNIT)2519:       CLOSE(LUNIT)
2699:    ENDIF2520:    ENDIF
2700: ELSE IF (CONFILET) THEN 2521: ELSE IF (CONFILET) THEN 
3004:           IF (CP==MUPDATE) CP=02825:           IF (CP==MUPDATE) CP=0
3005:      ENDDO2826:      ENDDO
3006:               2827:               
3007:      STP(1:D) = 1.0D02828:      STP(1:D) = 1.0D0
3008: ENDIF MAIN2829: ENDIF MAIN
3009: 2830: 
3010: !  Store the new search direction2831: !  Store the new search direction
3011: IF (NITERDONE.GT.1) SEARCHSTEP(POINT,1:D)=GTMP(1:D)2832: IF (NITERDONE.GT.1) SEARCHSTEP(POINT,1:D)=GTMP(1:D)
3012: 2833: 
3013: END SUBROUTINE MAKESTEP2834: END SUBROUTINE MAKESTEP
3014:  
3015: ! 
3016: ! Set up start and finish atom indices for each group once in intlbfgs 
3017: ! STARTGROUP(DOGROUP) and ENDGROUP(DOGROUP) 
3018: ! 
3019:  
3020: SUBROUTINE QCIOPTPERM(COORDSB,COORDSA,NATOMS,DEBUG,DOGROUP,PERM,STARTGROUP,IM) 
3021: USE KEY,ONLY : NPERMGROUP, NPERMSIZE, PERMGROUP, NSETS, SETS 
3022: IMPLICIT NONE 
3023: INTEGER PERM(NATOMS), DOGROUP, J1, J2, NATOMS, STARTGROUP(NPERMGROUP), J3, I1, I2, I3, IM 
3024: DOUBLE PRECISION COORDSA(3*NATOMS), COORDSB(3*NATOMS), D123, D132, D231, D213, D312, D321, DMIN, D12, D21 
3025: LOGICAL DEBUG, IDENTITY 
3026:  
3027: IDENTITY=.TRUE. 
3028: DO J1=1,NATOMS 
3029:    PERM(J1)=J1 
3030: ENDDO 
3031:  
3032: IF (NPERMSIZE(DOGROUP).EQ.3) THEN 
3033:    I1=PERMGROUP(STARTGROUP(DOGROUP)) 
3034:    I2=PERMGROUP(STARTGROUP(DOGROUP)+1) 
3035:    I3=PERMGROUP(STARTGROUP(DOGROUP)+2) 
3036: ! permutation 1 2 3 
3037:    D123=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I1-1)+3))**2 & 
3038:  &     +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I2-1)+3))**2 & 
3039:  &     +(COORDSB(3*(I3-1)+1)-COORDSA(3*(I3-1)+1))**2+(COORDSB(3*(I3-1)+2)-COORDSA(3*(I3-1)+2))**2+(COORDSB(3*(I3-1)+3)-COORDSA(3*(I3-1)+3))**2 
3040:    DMIN=D123 
3041: ! permutation 1 3 2 
3042:    D132=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I1-1)+3))**2 & 
3043:  &     +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I3-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I3-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I3-1)+3))**2 & 
3044:  &     +(COORDSB(3*(I3-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I3-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I3-1)+3)-COORDSA(3*(I2-1)+3))**2 
3045:    IF (D132.LT.DMIN) THEN  
3046:       DMIN=D132 
3047:       PERM(I1)=I1 
3048:       PERM(I2)=I3 
3049:       PERM(I3)=I2 
3050:       IDENTITY=.FALSE. 
3051:    ENDIF 
3052: ! permutation 2 1 3 
3053:    D213=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I2-1)+3))**2 & 
3054:  &     +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I1-1)+3))**2 & 
3055:  &     +(COORDSB(3*(I3-1)+1)-COORDSA(3*(I3-1)+1))**2+(COORDSB(3*(I3-1)+2)-COORDSA(3*(I3-1)+2))**2+(COORDSB(3*(I3-1)+3)-COORDSA(3*(I3-1)+3))**2 
3056:    IF (D213.LT.DMIN) THEN  
3057:       DMIN=D213 
3058:       PERM(I1)=I2 
3059:       PERM(I2)=I1 
3060:       PERM(I3)=I3 
3061:       IDENTITY=.FALSE. 
3062:    ENDIF 
3063: ! permutation 2 3 1 
3064:    D231=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I2-1)+3))**2 & 
3065:  &     +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I3-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I3-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I3-1)+3))**2 & 
3066:  &     +(COORDSB(3*(I3-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I3-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I3-1)+3)-COORDSA(3*(I1-1)+3))**2 
3067:    IF (D231.LT.DMIN) THEN  
3068:       DMIN=D231 
3069:       PERM(I1)=I2 
3070:       PERM(I2)=I3 
3071:       PERM(I3)=I1 
3072:       IDENTITY=.FALSE. 
3073:    ENDIF 
3074: ! permutation 3 2 1 
3075:    D321=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I3-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I3-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I3-1)+3))**2 & 
3076:  &     +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I2-1)+3))**2 & 
3077:  &     +(COORDSB(3*(I3-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I3-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I3-1)+3)-COORDSA(3*(I1-1)+3))**2 
3078:    IF (D321.LT.DMIN) THEN  
3079:       DMIN=D321 
3080:       PERM(I1)=I3 
3081:       PERM(I2)=I2 
3082:       PERM(I3)=I1 
3083:       IDENTITY=.FALSE. 
3084:    ENDIF 
3085: ! permutation 3 1 2 
3086:    D312=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I3-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I3-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I3-1)+3))**2 & 
3087:  &     +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I1-1)+3))**2 & 
3088:  &     +(COORDSB(3*(I3-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I3-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I3-1)+3)-COORDSA(3*(I2-1)+3))**2 
3089:    IF (D312.LT.DMIN) THEN  
3090:       DMIN=D312 
3091:       PERM(I1)=I3 
3092:       PERM(I2)=I1 
3093:       PERM(I3)=I2 
3094:       IDENTITY=.FALSE. 
3095:    ENDIF 
3096: !  IF (.NOT.IDENTITY) WRITE(*,'(A,2I6,6G15.5,3I6)') ' qcioptperm> images,D123,D132,D213,D231,D312,D321,permutation: ', & 
3097: !&               IM,IM+1,D123,D132,D213,D231,D312,D321,PERM(I1),PERM(I2),PERM(I3) 
3098:    WRITE(*,'(A,2I6,6G15.5,3I6)') ' qcioptperm> images,D123,D132,D213,D231,D312,D321,permutation: ', & 
3099:  &               IM,IM+1,D123,D132,D213,D231,D312,D321,PERM(I1),PERM(I2),PERM(I3) 
3100: ELSE IF (NPERMSIZE(DOGROUP).EQ.2) THEN 
3101: ! permutation 1 2 
3102:    I1=PERMGROUP(STARTGROUP(DOGROUP)) 
3103:    I2=PERMGROUP(STARTGROUP(DOGROUP)+1) 
3104:    D12=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I1-1)+3))**2 & 
3105:  &    +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I2-1)+3))**2  
3106:    DMIN=D12 
3107: ! permutation 2 1 
3108:    D21=(COORDSB(3*(I1-1)+1)-COORDSA(3*(I2-1)+1))**2+(COORDSB(3*(I1-1)+2)-COORDSA(3*(I2-1)+2))**2+(COORDSB(3*(I1-1)+3)-COORDSA(3*(I2-1)+3))**2 & 
3109:  &    +(COORDSB(3*(I2-1)+1)-COORDSA(3*(I1-1)+1))**2+(COORDSB(3*(I2-1)+2)-COORDSA(3*(I1-1)+2))**2+(COORDSB(3*(I2-1)+3)-COORDSA(3*(I1-1)+3))**2  
3110:    IF (D21.LT.DMIN) THEN  
3111:       DMIN=D21 
3112:       PERM(I1)=I2 
3113:       PERM(I2)=I1 
3114:       IF (NSETS(DOGROUP).GT.0) THEN 
3115:          DO J3=1,NSETS(DOGROUP) ! can be 2 or 3 
3116: !           WRITE(*,'(A,4I6)') 'I1,I2,sets(I1,j3),sets(I2,j3)=',I1,I2,sets(I1,j3),sets(I2,j3) 
3117:             PERM(SETS(I1,J3))=SETS(I2,J3) 
3118:             PERM(SETS(I2,J3))=SETS(I1,J3) 
3119:          ENDDO 
3120:       ENDIF 
3121:  
3122:    ENDIF 
3123:    IF (.NOT.IDENTITY) THEN 
3124:       WRITE(*,'(A,2I6,2G15.5,2I6)') ' qcioptperm> images,D12,D21,permutation: ',IM,IM+1,D12,D21,PERM(I1),PERM(I2) 
3125:       IF (NSETS(DOGROUP).GT.0) WRITE(*,'(A,6I6)') ' qcioptperm> associated permutation: ', & 
3126:   &                         (PERM(SETS(I1,J3)),PERM(SETS(I2,J3)),J3=1,NSETS(DOGROUP))   
3127:    ENDIF 
3128: ELSE 
3129:    WRITE(*,'(A,I6,A,I6)') ' qcioptperm> Unknown group size ',NPERMSIZE(DOGROUP),' for group ',DOGROUP 
3130:    STOP 
3131: ENDIF 
3132:  
3133: END SUBROUTINE QCIOPTPERM 


r33407/key.f90 2017-10-23 10:30:21.481312983 +0100 r33406/key.f90 2017-10-23 10:30:23.165335261 +0100
 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, & 19:      &        NCONSTRAINTFIX, INTIMAGECHECK, NREPULSIVEFIX, NRANROT, NENDDUP, LOCALPERMNEIGH, &
 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, & 20:      &        LOCALPERMMAXSEP, NONEDAPBC, STRUC, QCHEMESNAO, QCHEMESNMO, QCHEMESNZERO, QCHEMESNELEC, PMPATHINR, &
 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, & 21:      &        MULTISUNIT, MULTIFUNIT,NIMAGEINST,NGLJ,ST_TSSTEP,LANSTEP,NONFREEZE, &
 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, & 22:      &        MCPATHBINS,MCPATHEQUIL,MCPATHSTEPS,MCPATHPRTFRQ,MCPATHTS,MCPATHSCHECK,RPHSLICES,RPHQBINS, &
 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, & 23:      &        ITWIST, JTWIST, KTWIST, LTWIST, MCPATHSTART, MCPATHBLOCK, MCPATHOVER, NCPU, MCPATHDOBLOCK, MCMERGES, MCMERGEQ, &
 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, & 24:      &        MCMERGEI,GAUSSIANCHARGE,GAUSSIANMULTI,ITG03, REDOTS, QCIPERMCHECKINT, &
 25:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS, & 25:      &        MLPIN, MLPSTART, MLPOUT, MLPHIDDEN, MLPDATA, NMLP, N_TO_ALIGN, DJWRBID, STM, NHEXAMERS, &
 26:      &        MLQIN, MLQSTART, MLQOUT, MLQDATA, NMLQ, & 26:      &        MLQIN, MLQSTART, MLQOUT, MLQDATA, NMLQ, &
 27:      &        QCIADDREP, QCIBONDS, QCISECOND, MAXNACTIVE, QCIIMAGE, NADDTARGET, NUMNN, MULTI_COUNT, MULTI_LAST, MULTI_STEP, & 27:      &        QCIADDREP, QCIBONDS, QCISECOND, MAXNACTIVE, QCIIMAGE, NADDTARGET, NUMNN, MULTI_COUNT, MULTI_LAST, MULTI_STEP, &
 28:      &        NDOF, RECCOUNT, MLPPROBPOS, PUSHOPTMAX, MLPNEIGH, QCICYCN, QCIPDINT, NPEAKS, NDISPLACEMENTS, NROTATIONS, & 28:      &        NDOF, RECCOUNT, MLPPROBPOS, PUSHOPTMAX, MLPNEIGH, QCICYCN, QCIPDINT, NPEAKS, NDISPLACEMENTS, NROTATIONS, &
 29:      &        MAX_ANGMOM, BNB_NSTEPS, QUIPZ 29:      &        MAX_ANGMOM, BNB_NSTEPS
 30:  30: 
 31:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, & 31:       LOGICAL :: DTEST, MASST, RTEST, EFSTEPST, VECTORST, SUMMARYT, DUMPV, DUMPMAG, FREEZE, FREEZERANGE, GRADSQ, &
 32:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, & 32:      &        PGRAD, VALUEST, ADMT, BFGSMINT, BFGSTST, CHECKINDEX, TOSI, CONTAINER, &
 33:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, & 33:      &        GAUSSIAN, CADPAC, PRESSURE, FTEST, DCHECK, CP2K, DFTP, CPMD, CPMDC, FREEZERES, DF1T, &
 34:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, & 34:      &        VARIABLES, FIELDT, OHT, IHT, TDT, D5HT, TWOENDS, PV, FRACTIONAL, BLNT, HYBRIDMINT, &
 35:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, & 35:      &        INDEXT, LANCZOST, NOSHIFT, GAMESSUS, GAMESSUK, PVTS, RIGIDBODY, CASTEP, ONETEP, QCHEM, QCHEMES, VASP, &
 36:      &        BFGSSTEP, EFOLSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, & 36:      &        BFGSSTEP, EFOLSTEP, BULKT, HUPDATE, NOHESS, READV, NOIT, THOMSONT, SIO2T, SIO2C6T, BISECTT, BISECTDEBUG, &
 37:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, & 37:      &        TOSIC6, TOSIPOL, FIXIMAGE, DFTBT, CHECKCONT, CHECKDT, SHIFTED, READSP, DUMPSP, NOFRQS, &
 38:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, & 38:      &        ALLSTEPS, ALLVECTORS, MWVECTORS, WELCH, BINARY, READHESS, MOVIE, NORESET, TWOD, &
 39:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,& 39:      &        DOUBLET, REOPT, PARALLEL, LINEMIN, FIXD, KEEPINDEX, BSMIN, PRINTPTS, RKMIN, REPELTST,&
 53:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 53:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 54:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 54:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 55:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 55:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 56:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 56:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 57:      &        PBST, SSHT, GAUSSIAN03, GAUSSIAN09, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 57:      &        PBST, SSHT, GAUSSIAN03, GAUSSIAN09, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 58:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 58:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 59:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 59:      &        MALONALDEHYDE, SIO2PT, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 60:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 60:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 61:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T, & 61:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT, PYADD2T, SANDBOXT, LJADD3T, LJADD4T, &
 62:      &        MBPOLT, MULTIJOB_MACHINET, DUMPDATA_MACHINET, PLUSSIDET, MINUSSIDET, PUSHOPTT, MLPVB3NNT, GAUSSIAN16, QCICYCLEST, & 62:      &        MBPOLT, MULTIJOB_MACHINET, DUMPDATA_MACHINET, PLUSSIDET, MINUSSIDET, PUSHOPTT, MLPVB3NNT, GAUSSIAN16, QCICYCLEST, &
 63:      &        QCIDNEBT, QCIRESTART, QCILPERMDIST, FASTOVERLAPT, BNB_ALIGNT, QUIPT, QCIADDACIDT 63:      &        QCIDNEBT, QCIRESTART, QCILPERMDIST, FASTOVERLAPT, BNB_ALIGNT
 64:  64: 
 65:  65: 
 66: ! sy349 > for testing the flatpath after dneb 66: ! sy349 > for testing the flatpath after dneb
 67:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 67:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 68:       LOGICAL FLATPATHT 68:       LOGICAL FLATPATHT
 69:  69: 
 70: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 70: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 71:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 71:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 72:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 72:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 73:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 73:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
106:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &106:      &        REDOK, REDOFRAC, D1INIT, D2INIT, REDOE1, REDOE2, RPBETA, REPCON, PFORCE, &
107:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &107:      &        CPCONSTRAINTTOL, CPCONSTRAINTDEL, CPCONSTRAINTREP, CPCONSTRAINREPCUT, CPCONFRAC, &
108:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &108:      &        INTLJTOL, INTLJDEL, INTLJEPS, IMSEPMIN, IMSEPMAX, TRAPK, MINOVERLAP, &
109:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &109:      &        INTFREEZETOL, LOCALPERMCUT, LOCALPERMCUT2, LOCALPERMCUTINC, CHECKREPCUTOFF, CONCUTABS, &
110:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &110:      &        CONCUTFRAC, ENDNUMHESSDELTA, DNEBEFRAC, QCHEMSCALE, KAA, SIGMAAA, QUIPATOMMASS, TEMPERATURE1, &
111:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &111:      &        DISTORTINST,DELTAINST,MOLPROSCALE,COVER,STTSRMSCONV,LAN_DIST,LANCONV,LANFACTOR, &
112:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &112:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
113:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &113:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
114:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &114:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
115:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2, &115:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2, &
116:      &        TANHFAC, LJADDCUTOFF,LJADDREFNORM,MAXIMFACTOR, PUSHOPTCONV, QCICYCDIST, QCIPERMCUT, KERNELWIDTH, QUIPLATT(3,3)116:      &        TANHFAC, LJADDCUTOFF,LJADDREFNORM,MAXIMFACTOR, PUSHOPTCONV, QCICYCDIST, QCIPERMCUT, KERNELWIDTH
117: 117: 
118: !     sf344118: !     sf344
119:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &119:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
120:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &120:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
121:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT121:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT
122:  122:  
123:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)123:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:)
124:       DOUBLE PRECISION, ALLOCATABLE :: SHIFTL(:)124:       DOUBLE PRECISION, ALLOCATABLE :: SHIFTL(:)
125:       LOGICAL, ALLOCATABLE :: uniaxarray(:)125:       LOGICAL, ALLOCATABLE :: uniaxarray(:)
126:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)126:       DOUBLE PRECISION, ALLOCATABLE :: PYADDREP(:,:), PYADDATT(:,:)


r33407/keywords.f 2017-10-23 10:30:21.833317626 +0100 r33406/keywords.f 2017-10-23 10:30:23.393338277 +0100
644:          INTLJDEL=0.1D0644:          INTLJDEL=0.1D0
645:          INTLJEPS=1.0D0645:          INTLJEPS=1.0D0
646: 646: 
647: !647: !
648: ! QCI parameters648: ! QCI parameters
649: !649: !
650:          CONDATT=.FALSE.650:          CONDATT=.FALSE.
651:          QCIPOTT=.FALSE.651:          QCIPOTT=.FALSE.
652:          QCIPOT2T=.FALSE.652:          QCIPOT2T=.FALSE.
653:          QCIADDREP=0653:          QCIADDREP=0
654:          QCIADDACIDT=.FALSE. 
655:          QCIADDREPCUT=1.0D0654:          QCIADDREPCUT=1.0D0
656:          QCIADDREPEPS=1.0D0655:          QCIADDREPEPS=1.0D0
657:          QCINOREPINT=.FALSE.656:          QCINOREPINT=.FALSE.
658:          MAXNACTIVE=0657:          MAXNACTIVE=0
659: 658: 
660:          FREEZETOL=1.0D-3659:          FREEZETOL=1.0D-3
661:          FLATTESTT=.FALSE.660:          FLATTESTT=.FALSE.
662:          FLATEDIFF=1.0D-6661:          FLATEDIFF=1.0D-6
663:          QCIPERMCHECK=.FALSE.662:          QCIPERMCHECK=.FALSE.
664:          QCIPERMCHECKINT=100663:          QCIPERMCHECKINT=100
937:          DF1T=.FALSE.936:          DF1T=.FALSE.
938:          PULLT=.FALSE.937:          PULLT=.FALSE.
939:          CHEMSHIFT=.FALSE.938:          CHEMSHIFT=.FALSE.
940:          METRICTENSOR=.FALSE.939:          METRICTENSOR=.FALSE.
941: 940: 
942:          FRQCONV = 1.0D0941:          FRQCONV = 1.0D0
943:          DUMMY_FRQCONV = 0.0D0942:          DUMMY_FRQCONV = 0.0D0
944: 943: 
945:          QUIPARGSTRT=.FALSE.944:          QUIPARGSTRT=.FALSE.
946:          QUIPPARAMST=.FALSE.945:          QUIPPARAMST=.FALSE.
947:          QUIPZ=14 ! default is silicon 
948:          QUIPT=.FALSE. 
949: 946: 
950:          EIGENONLY=.FALSE.947:          EIGENONLY=.FALSE.
951:          COVER=0.990d0948:          COVER=0.990d0
952:          OVERCONV=.FALSE.949:          OVERCONV=.FALSE.
953:          TRUSTMODET=.FALSE.950:          TRUSTMODET=.FALSE.
954:          TMRATIO=0.71D0951:          TMRATIO=0.71D0
955: 952: 
956:          JPARAM=0.001D0953:          JPARAM=0.001D0
957:          RPHT=.FALSE.954:          RPHT=.FALSE.
958:          MCPATHT=.FALSE.955:          MCPATHT=.FALSE.
1018: 1015: 
1019:          MBPOLT= .FALSE.1016:          MBPOLT= .FALSE.
1020: ! OPEP stuff1017: ! OPEP stuff
1021:          OPEPT = .FALSE.1018:          OPEPT = .FALSE.
1022:          OPEP_RNAT = .FALSE.1019:          OPEP_RNAT = .FALSE.
1023: 1020: 
1024: ! cs675 ORBITALS stuff1021: ! cs675 ORBITALS stuff
1025:          ORBITALS = .FALSE.1022:          ORBITALS = .FALSE.
1026:          ORBVAREXPONENT = -11023:          ORBVAREXPONENT = -1
1027: 1024: 
 1025: ! Align routine variables
 1026:          FASTOVERLAPT = .FALSE.
 1027:          BNB_ALIGNT = .FALSE.
 1028:          KERNELWIDTH = -1       ! This is set to a sensible value in the ALIGN subroutine itself.
 1029:          NPEAKS = 1
 1030:          NDISPLACEMENTS = 1
 1031:          NROTATIONS = 1
 1032:          MAX_ANGMOM = 15
 1033:          BNB_NSTEPS = 2000
 1034: 
 1035:          !!!!! END OF VARIABLE INITIALISATION, START READING THE DATA FILE HERE !!!!!
 1036: 
1028:          CLSTRINGT=.FALSE.1037:          CLSTRINGT=.FALSE.
1029:          CLSTRINGTST=.FALSE.1038:          CLSTRINGTST=.FALSE.
1030:          IF (FILTH2.EQ.0) THEN1039:          IF (FILTH2.EQ.0) THEN
1031:             OPEN (5,FILE='odata',STATUS='OLD')1040:             OPEN (5,FILE='odata',STATUS='OLD')
1032:          ELSE1041:          ELSE
1033:             WRITE(OTEMP,*) FILTH21042:             WRITE(OTEMP,*) FILTH2
1034:             WRITE(OSTRING,'(A)') 'odata.' // TRIM(ADJUSTL(OTEMP))1043:             WRITE(OSTRING,'(A)') 'odata.' // TRIM(ADJUSTL(OTEMP))
1035:             OPEN (5,FILE=OSTRING,STATUS='OLD')1044:             OPEN (5,FILE=OSTRING,STATUS='OLD')
1036:          ENDIF1045:          ENDIF
1037:          DATA_UNIT=51046:          DATA_UNIT=5
1755: ! 1764: ! 
1756: ! BSMIN calculates a steepest-descent path using gradient only information1765: ! BSMIN calculates a steepest-descent path using gradient only information
1757: ! with convergence criterion GMAX for the RMS force and initial precision1766: ! with convergence criterion GMAX for the RMS force and initial precision
1758: ! EPS. The Bulirsch-Stoer algorithm is used.1767: ! EPS. The Bulirsch-Stoer algorithm is used.
1759: ! 1768: ! 
1760:          ELSE IF (WORD.EQ.'BSMIN') THEN1769:          ELSE IF (WORD.EQ.'BSMIN') THEN
1761:             BSMIN=.TRUE.1770:             BSMIN=.TRUE.
1762:             IF (NITEMS.GT.1) CALL READF(GMAX)1771:             IF (NITEMS.GT.1) CALL READF(GMAX)
1763:             IF (NITEMS.GT.2) CALL READF(EPS)1772:             IF (NITEMS.GT.2) CALL READF(EPS)
1764: 1773: 
 1774:          ELSE IF (WORD.EQ.'BRANCHNBOUND') THEN
 1775:             BNB_ALIGNT = .TRUE.
 1776:             IF (NITEMS.GT.1) CALL READI(BNB_NSTEPS) ! Default 2000
 1777: 
1765:          ELSE IF (WORD.EQ.'BULK') THEN1778:          ELSE IF (WORD.EQ.'BULK') THEN
1766:             BULKT=.TRUE.1779:             BULKT=.TRUE.
1767:             IF (NITEMS.GT.1) THEN1780:             IF (NITEMS.GT.1) THEN
1768:                CALL READF(bulk_boxvec(1))1781:                CALL READF(bulk_boxvec(1))
1769:             ENDIF1782:             ENDIF
1770:             IF (NITEMS.GT.2) THEN1783:             IF (NITEMS.GT.2) THEN
1771:                CALL READF(bulk_boxvec(2))1784:                CALL READF(bulk_boxvec(2))
1772:             ENDIF1785:             ENDIF
1773:             IF (NITEMS.GT.3) THEN1786:             IF (NITEMS.GT.3) THEN
1774:                CALL READF(bulk_boxvec(3))1787:                CALL READF(bulk_boxvec(3))
2968: ! FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF2981: ! FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
2969: ! 2982: ! 
2970: ! 2983: ! 
2971: ! Distance dependent dielectric for Paul Mortenson`s amber2984: ! Distance dependent dielectric for Paul Mortenson`s amber
2972: ! Obsolete keyword now commented out to avoid accidental use2985: ! Obsolete keyword now commented out to avoid accidental use
2973: ! 2986: ! 
2974: ! ELSE IF (WORD.EQ.'FAKEWATER') THEN2987: ! ELSE IF (WORD.EQ.'FAKEWATER') THEN
2975: ! FAKEWATER=.TRUE.2988: ! FAKEWATER=.TRUE.
2976: ! WRITE (*,'(A)') ' SETTINGS Distance dependent dielectric will be used'2989: ! WRITE (*,'(A)') ' SETTINGS Distance dependent dielectric will be used'
2977: ! 2990: ! 
 2991: 
 2992:          ELSE IF (WORD.EQ.'FASTOVERLAP') THEN
 2993:             FASTOVERLAPT = .TRUE.
 2994:             
 2995:             IF (NITEMS .GT. 1) THEN ! Read a parameter to determine how many displacements/rotations to consider for bulk/cluster
 2996:                                  ! alignments respectively
 2997: 
 2998:                CALL READI(NPEAKS) ! NPEAKS doesn't get used any time other than here.
 2999:                NDISPLACEMENTS=NPEAKS ! Default = 1. This will be used if we have a periodic system
 3000:                NROTATIONS=NPEAKS ! Default = 1. This will be used if we have a cluster (but we don't know BULKT yet)
 3001:                
 3002:                IF (NITEMS .GT. 2) THEN
 3003:                   CALL READF(KERNELWIDTH) ! Default = -1 (which leads to a sensible value being chosen when the subroutine is called)
 3004:             
 3005:                   IF (NITEMS .GT.3) THEN
 3006:                      CALL READI(MAX_ANGMOM) ! Default = 15 (but larger values might be better for large systems)
 3007:                      IF(BULKT) WRITE(*,*) "keywords> Warning: maximum angular moment value explicitly specified but",
 3008:      &                                    "this parameter is not used for periodic systems."
 3009:                   ENDIF
 3010:                ENDIF
 3011:             ENDIF
 3012: 
 3013: 
2978: ! Integer variable to distinguish output files from parallel maiden jobs3014: ! Integer variable to distinguish output files from parallel maiden jobs
2979: ! 3015: ! 
2980:          ELSE IF (WORD.EQ.'FILTH') THEN3016:          ELSE IF (WORD.EQ.'FILTH') THEN
2981:             IF (FILTH.EQ.0) THEN3017:             IF (FILTH.EQ.0) THEN
2982:                CALL READI(FILTH)3018:                CALL READI(FILTH)
2983:             ELSE3019:             ELSE
2984:                WRITE(*,'(A)') 'WARNING **** FILTH keyword in odata was overridden by command line argument'3020:                WRITE(*,'(A)') 'WARNING **** FILTH keyword in odata was overridden by command line argument'
2985:             ENDIF3021:             ENDIF
2986: ! 3022: ! 
2987: ! Specifies that FIXIMAGE should be set permanently after step3023: ! Specifies that FIXIMAGE should be set permanently after step
3561: ! Maximum distance for constrained atoms3597: ! Maximum distance for constrained atoms
3562: !3598: !
3563:          ELSE IF ((WORD.EQ.'INTCONCUT').OR.(WORD.EQ.'QCICONCUT')) THEN3599:          ELSE IF ((WORD.EQ.'INTCONCUT').OR.(WORD.EQ.'QCICONCUT')) THEN
3564:             CALL READF(INTCONCUT)3600:             CALL READF(INTCONCUT)
3565: !3601: !
3566: ! QCI LOPERMDIST checks for images3602: ! QCI LOPERMDIST checks for images
3567: !3603: !
3568:          ELSE IF ((WORD.EQ.'INTLPERMDIST').OR.(WORD.EQ.'QCILPERMDIST')) THEN3604:          ELSE IF ((WORD.EQ.'INTLPERMDIST').OR.(WORD.EQ.'QCILPERMDIST')) THEN
3569:             QCILPERMDIST=.TRUE.3605:             QCILPERMDIST=.TRUE.
3570:             IF (NITEMS.GT.1) CALL READI(QCIPDINT)   ! interval for check3606:             IF (NITEMS.GT.1) CALL READI(QCIPDINT)   ! interval for check
3571:             IF (NITEMS.GT.2) CALL READF(QCIPERMCUT) ! maximum allowed distance for valid alignment3607:             IF (NITEMS.GT.2) CALL READI(QCIPERMCUT) ! maximum allowed distance for valid alignment
3572:             WRITE(*,'(A,I8,A,G20.10)') ' keywords> QCI active permutational alignment check every ',QCIPDINT,' steps, tolerance=', 3608:             WRITE(*,'(A,I8,A,G20.10)') ' keywords> QCI active permutational alignment check every ',QCIPDINT,' steps, tolerance=', 
3573:      &                                  QCIPERMCUT     3609:      &                                  QCIPERMCUT     
3574: ! 3610: ! 
3575: ! Use constraint potential for initial interpolation in each cycle.3611: ! Use constraint potential for initial interpolation in each cycle.
3576: ! 3612: ! 
3577:          ELSE IF ((WORD.EQ.'INTCONSTRAINT').OR.(WORD.EQ.'QCI')) THEN3613:          ELSE IF ((WORD.EQ.'INTCONSTRAINT').OR.(WORD.EQ.'QCI')) THEN
3578:             INTCONSTRAINTT=.TRUE.3614:             INTCONSTRAINTT=.TRUE.
3579:             IF (NITEMS.GT.1) CALL READF(INTCONSTRAINTTOL)3615:             IF (NITEMS.GT.1) CALL READF(INTCONSTRAINTTOL)
3580:             IF (NITEMS.GT.2) CALL READF(INTCONSTRAINTDEL)3616:             IF (NITEMS.GT.2) CALL READF(INTCONSTRAINTDEL)
3581:             IF (NITEMS.GT.3) CALL READF(INTCONSTRAINTREP)3617:             IF (NITEMS.GT.3) CALL READF(INTCONSTRAINTREP)
3646:                   CLOSE(LUNIT)3682:                   CLOSE(LUNIT)
3647:                   PRINT '(A,I6,A)',' keyword> Read ',NCONGEOM,' reference geometries from congeom file'3683:                   PRINT '(A,I6,A)',' keyword> Read ',NCONGEOM,' reference geometries from congeom file'
3648:                   IF (NCONGEOM.LT.2) PRINT '(A)',' WARNING *** insufficient reference geometries - using end point minima'3684:                   IF (NCONGEOM.LT.2) PRINT '(A)',' WARNING *** insufficient reference geometries - using end point minima'
3649:                ENDIF3685:                ENDIF
3650:             ENDIF3686:             ENDIF
3651: ! 3687: ! 
3652: ! DO NOT Use the quasi-continuous metric for connection attempts, instead of distance.3688: ! DO NOT Use the quasi-continuous metric for connection attempts, instead of distance.
3653: ! 3689: ! 
3654:             INTERPCOSTFUNCTION=.FALSE.3690:             INTERPCOSTFUNCTION=.FALSE.
3655: !3691: !
3656: ! Add complete amino acids 
3657: ! 
3658:       ELSE IF (WORD.EQ.'QCIADDACID') THEN 
3659:          QCIADDACIDT=.TRUE. 
3660: ! 
3661: !Use topology information for QCI constraints for AMBER3692: !Use topology information for QCI constraints for AMBER
3662: !3693: !
3663:       ELSE IF (WORD.EQ.'QCIAMBER') THEN3694:       ELSE IF (WORD.EQ.'QCIAMBER') THEN
3664:          QCIAMBERT=.TRUE.3695:          QCIAMBERT=.TRUE.
3665:          WRITE(*,'(A)') ' keyword> Use topology file for constraints in QCI'3696:          WRITE(*,'(A)') ' keyword> Use topology file for constraints in QCI'
3666: 3697: 
3667:          ELSE IF (WORD.EQ.'INTFREEZE') THEN3698:          ELSE IF (WORD.EQ.'INTFREEZE') THEN
3668:             INTFREEZET=.TRUE.3699:             INTFREEZET=.TRUE.
3669:             IF (NITEMS.GT.1) CALL READF(INTFREEZETOL)3700:             IF (NITEMS.GT.1) CALL READF(INTFREEZETOL)
3670:             IF (NITEMS.GT.2) CALL READI(INTFREEZEMIN)3701:             IF (NITEMS.GT.2) CALL READI(INTFREEZEMIN)
6092: ! Coded by Javier.6123: ! Coded by Javier.
6093: ! 6124: ! 
6094:          ELSE IF (WORD.EQ.'QTIP4PF') THEN6125:          ELSE IF (WORD.EQ.'QTIP4PF') THEN
6095:             QTIP4PFT=.TRUE.6126:             QTIP4PFT=.TRUE.
6096: 6127: 
6097:             ! sn402: added (see comments at keyword FRQCONV)6128:             ! sn402: added (see comments at keyword FRQCONV)
6098:             FRQCONV = 2.045483D136129:             FRQCONV = 2.045483D13
6099:             WRITE(*,*) "keywords> Set frequency conversion factor to the QTIP4PF default value: ", FRQCONV6130:             WRITE(*,*) "keywords> Set frequency conversion factor to the QTIP4PF default value: ", FRQCONV
6100:             WRITE(*,*) "keywords> This corresponds to frequencies being given in radians/s"6131:             WRITE(*,*) "keywords> This corresponds to frequencies being given in radians/s"
6101: 6132: 
6102:          ELSE IF (WORD.EQ.'QUIPZ') THEN6133:          ELSE IF (WORD.EQ.'QUIPARGSTR') THEN
6103:             CALL READI(QUIPZ)6134:             QUIPARGSTRT=.TRUE.
 6135:             QARGSTR='IP LJ'
 6136:             CALL READA(QARGSTR)
 6137: 
 6138:          ELSE IF (WORD.EQ.'QUIPPARAMS') THEN
 6139:             QUIPPARAMST=.TRUE.
 6140:             QUIPATOMTYPE='Ag '
 6141:             QUIPATOMMASS=1.0D0
 6142:             IF (NITEMS.GT.1) THEN
 6143:                CALL READA(QUIPATOMTYPE)
 6144:                CALL READF(QUIPATOMMASS)
 6145:             ENDIF
6104: 6146: 
6105:          ELSE IF (WORD.EQ.'QUIPLATT') THEN 
6106:              QUIPT=.TRUE. 
6107:              CALL READF(QUIPLATT(1,1)) 
6108:              CALL READF(QUIPLATT(2,1)) 
6109:              CALL READF(QUIPLATT(3,1)) 
6110:              CALL READF(QUIPLATT(1,2)) 
6111:              CALL READF(QUIPLATT(2,2)) 
6112:              CALL READF(QUIPLATT(3,2)) 
6113:              CALL READF(QUIPLATT(1,3)) 
6114:              CALL READF(QUIPLATT(2,3)) 
6115:              CALL READF(QUIPLATT(3,3)) 
6116: ! 6147: ! 
6117: ! RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR6148: ! RRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR
6118: ! 6149: ! 
6119: ! 6150: ! 
6120: ! Spherical container6151: ! Spherical container
6121: ! 6152: ! 
6122:          ELSE IF (WORD.EQ.'RADIUS') THEN6153:          ELSE IF (WORD.EQ.'RADIUS') THEN
6123:             CONTAINER=.TRUE.6154:             CONTAINER=.TRUE.
6124:             CALL READF(RADIUS)6155:             CALL READF(RADIUS)
6125:             RADIUS=RADIUS**26156:             RADIUS=RADIUS**2


r33407/lopermdist.f90 2017-10-23 10:30:22.053320550 +0100 r33406/lopermdist.f90 2017-10-23 10:30:23.613341187 +0100
165: ! DMEAN, SORTLIST, TRIED, PERMUTABLE, and DLIST entries refer to original165: ! DMEAN, SORTLIST, TRIED, PERMUTABLE, and DLIST entries refer to original
166: ! atom labels. Use NEWPERM to find where they are in coordinate lists.166: ! atom labels. Use NEWPERM to find where they are in coordinate lists.
167: !167: !
168:    outer1: DO J2=1,NATOMS168:    outer1: DO J2=1,NATOMS
169: !169: !
170: ! Don't allow members of the same permutational group 170: ! Don't allow members of the same permutational group 
171: ! to appear as reference neighbours.171: ! to appear as reference neighbours.
172: !172: !
173:       IF (TRIED(J2).EQ.-1) THEN173:       IF (TRIED(J2).EQ.-1) THEN
174:          XDUMMY=1.0D9174:          XDUMMY=1.0D9
 175:       ELSE IF ((DOGROUP.GT.0).AND.(.NOT.ATOMACTIVE(J2))) THEN ! only use active atoms for QCI single group
 176:          XDUMMY=1.0D9
175:       ELSE177:       ELSE
176:          IF ((DOGROUP.GT.0).AND.(.NOT.ATOMACTIVE(J2))) THEN ! only use active atoms for QCI single group178:          DA=(XA-DUMMYA(3*(NEWPERM(J2)-1)+1))**2 &
177:             XDUMMY=1.0D9179:   &        +(YA-DUMMYA(3*(NEWPERM(J2)-1)+2))**2 &
178:          ELSE180:   &        +(ZA-DUMMYA(3*(NEWPERM(J2)-1)+3))**2
179:             DA=(XA-DUMMYA(3*(NEWPERM(J2)-1)+1))**2 &181: !        DB=(XB-DUMMYB(3*(NEWPERM(J2)-1)+1))**2 &
180:   &           +(YA-DUMMYA(3*(NEWPERM(J2)-1)+2))**2 &182: ! &        +(YB-DUMMYB(3*(NEWPERM(J2)-1)+2))**2 &
181:   &           +(ZA-DUMMYA(3*(NEWPERM(J2)-1)+3))**2183: ! &        +(ZB-DUMMYB(3*(NEWPERM(J2)-1)+3))**2
182: !           DB=(XB-DUMMYB(3*(NEWPERM(J2)-1)+1))**2 &184:          DB=(XB-DUMMYB(3*(J2-1)+1))**2 &
183: ! &           +(YB-DUMMYB(3*(NEWPERM(J2)-1)+2))**2 &185:   &        +(YB-DUMMYB(3*(J2-1)+2))**2 &
184: ! &           +(ZB-DUMMYB(3*(NEWPERM(J2)-1)+3))**2186:   &        +(ZB-DUMMYB(3*(J2-1)+3))**2
185:             DB=(XB-DUMMYB(3*(J2-1)+1))**2 &187:          XDUMMY=(SQRT(DA)+SQRT(DB))/2.0D0
186:   &           +(YB-DUMMYB(3*(J2-1)+2))**2 & 
187:   &           +(ZB-DUMMYB(3*(J2-1)+3))**2 
188:             XDUMMY=(SQRT(DA)+SQRT(DB))/2.0D0 
189:          ENDIF 
190:       ENDIF188:       ENDIF
191:       loop1: DO J3=1,J2189:       loop1: DO J3=1,J2
192:          IF (XDUMMY.LT.DMEAN(J3)) THEN190:          IF (XDUMMY.LT.DMEAN(J3)) THEN
193: !191: !
194: ! Move the rest down.192: ! Move the rest down.
195: !193: !
196:             DO J4=J2,J3+1,-1194:             DO J4=J2,J3+1,-1
197:                DMEAN(J4)=DMEAN(J4-1)195:                DMEAN(J4)=DMEAN(J4-1)
198:                SORTLIST(J4)=SORTLIST(J4-1)196:                SORTLIST(J4)=SORTLIST(J4-1)
199:             ENDDO197:             ENDDO
211:    LDBESTATOM=1.0D100209:    LDBESTATOM=1.0D100
212:    NOTHER=0210:    NOTHER=0
213:    DO J2=1,NATOMS211:    DO J2=1,NATOMS
214:       IF (TRIED(J2).EQ.1) THEN212:       IF (TRIED(J2).EQ.1) THEN
215:          NOTHER=NOTHER+1213:          NOTHER=NOTHER+1
216:          DLIST(NOTHER)=J2214:          DLIST(NOTHER)=J2
217:       ENDIF215:       ENDIF
218:    ENDDO216:    ENDDO
219:    ADDED=.FALSE.217:    ADDED=.FALSE.
220:    outer2: DO J2=1,NATOMS218:    outer2: DO J2=1,NATOMS
221: !     WRITE(*,'(A,I6,2G20.10)') 'lopermdist> J2,dmean,localpermcut2=',J2,DMEAN(J2),LOCALPERMCUT2 
222:       IF (DMEAN(J2).GT.LOCALPERMCUT2) THEN219:       IF (DMEAN(J2).GT.LOCALPERMCUT2) THEN
223: !        PRINT '(A)',' lopermdist> No more atoms within cutoff'220: !        PRINT '(A)',' lopermdist> No more atoms within cutoff'
224:          GOTO 91221:          GOTO 91
225:       ENDIF222:       ENDIF
226:       IF (TRIED(SORTLIST(J2)).EQ.0) THEN223:       IF (TRIED(SORTLIST(J2)).EQ.0) THEN
227:          ADDED=.TRUE.224:          ADDED=.TRUE.
228:          NOTHER=NOTHER+1225:          NOTHER=NOTHER+1
229:          IF (NOTHER+PATOMS.GT.NATOMS) THEN226:          IF (NOTHER+PATOMS.GT.NATOMS) THEN
230:             PRINT '(A,I6)', &227:             PRINT '(A,I6)', &
231:   & ' lopermdist> ERROR *** number of neighbours plus number of permutable atoms exceeds total for group ',J1228:   & ' lopermdist> ERROR *** number of neighbours plus number of permutable atoms exceeds total for group ',J1
362: !  PRINT '(A,4I6,2G20.10)','NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,root LDISTANCE,root LDBEST=', &359: !  PRINT '(A,4I6,2G20.10)','NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,root LDISTANCE,root LDBEST=', &
363: ! &                         NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,SQRT(LDISTANCE),SQRT(LDBEST(J1))360: ! &                         NORBIT1,NORBIT2,NCHOOSE1,NCHOOSE2,SQRT(LDISTANCE),SQRT(LDBEST(J1))
364: 361: 
365:    IF (NCHOOSE2.LT.NORBIT2) GOTO 30362:    IF (NCHOOSE2.LT.NORBIT2) GOTO 30
366:    IF (NCHOOSE1.LT.NORBIT1) GOTO 65363:    IF (NCHOOSE1.LT.NORBIT1) GOTO 65
367:    IF (NCHOOSEB2.LT.NORBITB2) GOTO 31364:    IF (NCHOOSEB2.LT.NORBITB2) GOTO 31
368:    IF (NCHOOSEB1.LT.NORBITB1) GOTO 66365:    IF (NCHOOSEB1.LT.NORBITB1) GOTO 66
369: 366: 
370: !  PRINT '(A,2G20.10)','root LDBESTATOM,LOCALPERMCUT=',SQRT(LDBESTATOM),LOCALPERMCUT 367: !  PRINT '(A,2G20.10)','root LDBESTATOM,LOCALPERMCUT=',SQRT(LDBESTATOM),LOCALPERMCUT 
371:    IF (SQRT(LDBESTATOM).GT.LOCALPERMCUT) THEN368:    IF (SQRT(LDBESTATOM).GT.LOCALPERMCUT) THEN
372: !     IF (DEBUG) THEN369:       IF (DEBUG) THEN
373: !        PRINT '(A,G15.5,A,I6,A,G15.5)',' lopermdist> Best distance ',SQRT(LDBESTATOM), &370: !        PRINT '(A,G15.5,A,I6)',' lopermdist> Best distance ',SQRT(LDBESTATOM), &
374: ! &                                     ' is too large for atom ',DLIST(NOTHER),' cutoff=',LOCALPERMCUT371: ! &                                     ' is too large for atom ',DLIST(NOTHER)
375: !     ENDIF372:       ENDIF
376:       TRIED(DLIST(NOTHER))=-1373:       TRIED(DLIST(NOTHER))=-1
377:       IF (NADDED.GT.1) THEN374:       IF (NADDED.GT.1) THEN
378: !        IF (DEBUG) THEN375:          IF (DEBUG) THEN
379: !           PRINT '(A)',' lopermdist> and partner atoms:'376: !           PRINT '(A)',' lopermdist> and partner atoms:'
380: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)377: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)
381: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:'378: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:'
382: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1)379: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1)
383: !        ENDIF380:          ENDIF
384:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1381:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=-1
385:       ENDIF382:       ENDIF
386:       GOTO 71383:       GOTO 71
387:    ELSE384:    ELSE
388: !     IF (DEBUG) PRINT '(A,G20.10,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), &385:       IF (DEBUG) PRINT '(A,G20.10,3(A,I6))',' lopermdist> Best distance ',SQRT(LDBESTATOM), &
389: ! &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours' 386:   &                    ' is OK for myorient with atom ',DLIST(NOTHER),' and ',NOTHER,' neighbours' 
390:       TRIED(DLIST(NOTHER))=1387:       TRIED(DLIST(NOTHER))=1
391:       IF (NADDED.GT.1) THEN388:       IF (NADDED.GT.1) THEN
392: !        IF (DEBUG) THEN389: !        IF (DEBUG) THEN
393: !           PRINT '(A)',' lopermdist> and partner atoms:'390: !           PRINT '(A)',' lopermdist> and partner atoms:'
394: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)391: !           PRINT '(20I5)',DLIST(NOTHER-NADDED+1:NOTHER-1)
395: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:'392: !           IF (DOGROUP.GT.0) PRINT '(A)',' lopermdist> atomactive values:'
396: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1)393: !           IF (DOGROUP.GT.0) PRINT '(20L5)',ATOMACTIVE(NOTHER-NADDED+1:NOTHER-1)
397: !        ENDIF394: !        ENDIF
398:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1395:          TRIED(DLIST(NOTHER-NADDED+1:NOTHER-1))=1
399:       ENDIF396:       ENDIF
447:          DO J3=1,NSETS(J1)444:          DO J3=1,NSETS(J1)
448:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3)445:             SAVEPERM(SETS(PERMGROUP(NDUMMY+J2-1),J3))=SETS(NEWPERM(PERMGROUP(NDUMMY+LPERM(J2)-1)),J3)
449:          ENDDO446:          ENDDO
450:       ENDDO447:       ENDDO
451:    ENDIF448:    ENDIF
452: !449: !
453: ! Save current optimal permutation in NEWPERM450: ! Save current optimal permutation in NEWPERM
454: !451: !
455:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS)452:    NEWPERM(1:NATOMS)=SAVEPERM(1:NATOMS)
456:    DSUM=DSUM+SQRT(LDBEST(J1))453:    DSUM=DSUM+SQRT(LDBEST(J1))
457: !  PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' best distance=',SQRT(LDBEST(J1)),' total=',DSUM454:    PRINT '(A,I6,2(A,F20.10))',' lopermdist> For group ',J1,' best distance=',SQRT(LDBEST(J1)),' total=',DSUM
458: !  PRINT '(A)','best permutation is now'455: !  PRINT '(A)','best permutation is now'
459: !  PRINT '(20I6)',NEWPERM(1:NATOMS)456: !  PRINT '(20I6)',NEWPERM(1:NATOMS)
460: 457: 
461: !458: !
462: ! Update NDUMMY, the cumulative offset for PERMGROUP459: ! Update NDUMMY, the cumulative offset for PERMGROUP
463: !460: !
464: 864   NDUMMY=NDUMMY+NPERMSIZE(J1)461: 864   NDUMMY=NDUMMY+NPERMSIZE(J1)
465: ENDDO  !  end of loop over groups of permutable atoms462: ENDDO  !  end of loop over groups of permutable atoms
466: 463: 
467: IF (DOGROUP.GT.0) THEN464: IF (DOGROUP.GT.0) THEN
468:    DISTANCE=SQRT(LDBEST(DOGROUP))465:    DISTANCE=SQRT(LDBEST(DOGROUP))
469:    NMOVE=0466:    NMOVE=0
470:    DO J2=1,NATOMS467:    DO J2=1,NATOMS
471:       IF (NEWPERM(J2).NE.J2) THEN468:       IF (NEWPERM(J2).NE.J2) THEN
472: !        WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2469: !        WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2
473:          NMOVE=NMOVE+1470:          NMOVE=NMOVE+1
474:       ENDIF471:       ENDIF
475:    ENDDO472:    ENDDO
476: !  PRINT '(A,I6,A,I6)',' lopermdist> Total permutations (not applied!) for optimal alignment of group ',DOGROUP,' is ',NMOVE473: !  PRINT '(A,I6,A,I6)',' lopermdist> Total permutations for optimal alignment of group ',DOGROUP,' is ',NMOVE
477:    RETURN474:    RETURN
478: ENDIF475: ENDIF
479: NMOVE=0476: NMOVE=0
480: ! DO J2=1,NATOMS477: DO J2=1,NATOMS
481: !    IF (NEWPERM(J2).NE.J2) THEN478:    IF (NEWPERM(J2).NE.J2) THEN
482: !       WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2479: !     WRITE(*,'(A,I6,A,I6)') ' lopermdist> need to move atom ',NEWPERM(J2),' to position ',J2
483: !       NMOVE=NMOVE+1480:       NMOVE=NMOVE+1
484: !    ENDIF481:    ENDIF
485: ! ENDDO482: ENDDO
486:   IF (DEBUG) PRINT '(A,I6)',' lopermdist> Total permutations for optimal alignment (will be applied)=',NMOVE483: ! PRINT '(A,I6)',' lopermdist> Total permutations for optimal alignment=',NMOVE
487: !484: !
488: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA485: ! NEWPERM(J1) is the atom that moves to position J1 to map COORDSA
489: ! to the current best alignment. 486: ! to the current best alignment. 
490: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current487: ! This loop just appears to set SAVEPERM and ALLPERM equal to the current
491: ! NEWPERM.488: ! NEWPERM.
492: !489: !
493: !490: !
494: ! Putting the ALLPERM(J1)=J1 into the second loop causes pgf90 to miscompile!!491: ! Putting the ALLPERM(J1)=J1 into the second loop causes pgf90 to miscompile!!
495: !492: !
496: DO J1=1,NATOMS493: DO J1=1,NATOMS


r33407/make_conpot.f90 2017-10-23 10:30:22.273323460 +0100 r33406/make_conpot.f90 2017-10-23 10:30:23.833344096 +0100
 41: ! 41: !
 42: ! If this is not the first call, and we are being passed two minima, 42: ! If this is not the first call, and we are being passed two minima,
 43: ! then we are doing an interpolation metric for a new pair of minima. 43: ! then we are doing an interpolation metric for a new pair of minima.
 44: ! We should optimise the permutational isomers on reference minimum 1 44: ! We should optimise the permutational isomers on reference minimum 1
 45: ! and then do the overall alignment with newmindist, fixing the 45: ! and then do the overall alignment with newmindist, fixing the
 46: ! permutational isomers. This should put the permutational isomers 46: ! permutational isomers. This should put the permutational isomers
 47: ! in register with the constraints, which were calculated for all 47: ! in register with the constraints, which were calculated for all
 48: ! the reference minima after aligning with the first. 48: ! the reference minima after aligning with the first.
 49: ! 49: !
 50:    IF ((CALLED.OR.CONDATT).AND.(NCPFIT.EQ.2)) THEN 50:    IF ((CALLED.OR.CONDATT).AND.(NCPFIT.EQ.2)) THEN
 51:       CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),MINCOORDS(1,1:3*NATOMS),NATOMS,DEBUG, & 51:       CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),MINCOORDS(1,1:3*NATOMS),NATOMS,DEBUG, &
 52:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT) 52:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT)
 53:       CALL MINPERMDIST(CONGEOM(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DEBUG, & 53:       CALL ALIGN_DECIDE(CONGEOM(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DEBUG, &
 54:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT) 54:   &                    PARAM1,PARAM2,PARAM3,BULKT,TWOD,D,D2,RIGIDBODY,RMAT)
 55:       CALL NEWMINDIST(MINCOORDS(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DISTANCE, & 55:       CALL NEWMINDIST(MINCOORDS(1,1:3*NATOMS),MINCOORDS(2,1:3*NATOMS),NATOMS,DISTANCE, &
 56:   &                   BULKT,TWOD,'AX    ',.FALSE.,RIGIDBODY,DEBUG,RMAT) 56:   &                   BULKT,TWOD,'AX    ',.FALSE.,RIGIDBODY,DEBUG,RMAT)
 57:    ENDIF 57:    ENDIF
 58: ENDIF 58: ENDIF
 59:  59: 
 60: NQCIFREEZE=0 60: NQCIFREEZE=0
 61: IF (FREEZE) THEN 61: IF (FREEZE) THEN
 62:    WRITE(*, '(A)') ' make_conpot> ERROR *** QCI has not been coded for frozen atoms yet' 62:    WRITE(*, '(A)') ' make_conpot> ERROR *** QCI has not been coded for frozen atoms yet'
 63:    STOP 63:    STOP
256:       256:       
257:       NREPULSIVEFIX=NREPULSIVEFIX+1257:       NREPULSIVEFIX=NREPULSIVEFIX+1
258:       IF (NREPULSIVEFIX.GT.NREPMAX) CALL REPDOUBLE258:       IF (NREPULSIVEFIX.GT.NREPMAX) CALL REPDOUBLE
259:       IF (NCONGEOM.GE.2) THEN259:       IF (NCONGEOM.GE.2) THEN
260:          REPIFIX(NREPULSIVEFIX)=J1260:          REPIFIX(NREPULSIVEFIX)=J1
261:          REPJFIX(NREPULSIVEFIX)=J2261:          REPJFIX(NREPULSIVEFIX)=J2
262:          REPCUTFIX(NREPULSIVEFIX)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)262:          REPCUTFIX(NREPULSIVEFIX)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)
263:       ENDIF263:       ENDIF
264:       REPI(NREPULSIVEFIX)=J1264:       REPI(NREPULSIVEFIX)=J1
265:       REPJ(NREPULSIVEFIX)=J2265:       REPJ(NREPULSIVEFIX)=J2
266:  
267: !266: !
268: ! Use the minimum of the end point distances and INTCONSTRAINREPCUT for each contact.267: ! Use the minimum of the end point distances and INTCONSTRAINREPCUT for each contact.
269: !268: !
270:       REPCUT(NREPULSIVEFIX)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)269:       REPCUT(NREPULSIVEFIX)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)
271:       IF (DEBUG) WRITE(*, '(A,I6,A,I6,A,F15.5,A,I10)') ' make_conpot> Adding repulsion for atom ',J1, &270:       IF (DEBUG) WRITE(*, '(A,I6,A,I6,A,F15.5,A,I10)') ' make_conpot> Adding repulsion for atom ',J1, &
272:   &              ' with atom ',J2,' cutoff=',DMIN,' # repulsions ',NREPULSIVEFIX271:   &              ' with atom ',J2,' cutoff=',DMIN,' # repulsions ',NREPULSIVEFIX
273: !     IF ((REPIFIX(NREPULSIVEFIX).EQ.2024).OR.(REPJFIX(NREPULSIVEFIX).EQ.2024)) THEN 
274: !        WRITE(*,'(A,I6,I6,3G20.10)') 'make_conpot> NREPULSIVEFIX,repifix,repjfix,df,dmin,REPCUTFIX=', & 
275: ! &         NREPULSIVEFIX,repifix(NREPULSIVEFIX),repjfix(NREPULSIVEFIX),df,dmin,REPCUTFIX(NREPULSIVEFIX) 
276: !     ENDIF 
277:    ENDDO rep2272:    ENDDO rep2
278: ENDDO273: ENDDO
279: 963 CONTINUE274: 963 CONTINUE
280: NREPULSIVE=NREPULSIVEFIX275: NREPULSIVE=NREPULSIVEFIX
281: 276: 
282: WRITE(*, '(A,2I10,A,G20.10)') ' make_conpot> Total number of constraints and repulsions=', &277: WRITE(*, '(A,2I10,A,G20.10)') ' make_conpot> Total number of constraints and repulsions=', &
283:   &   NCONSTRAINT,NREPULSIVE,' for tolerance parameter ',LINTCONSTRAINTTOL278:   &   NCONSTRAINT,NREPULSIVE,' for tolerance parameter ',LINTCONSTRAINTTOL
284: 279: 
285: IF (ALLOCATED(CONACTIVE)) DEALLOCATE(CONACTIVE)280: IF (ALLOCATED(CONACTIVE)) DEALLOCATE(CONACTIVE)
286: ALLOCATE(CONACTIVE(NCONSTRAINT))281: ALLOCATE(CONACTIVE(NCONSTRAINT))


r33407/minpermdist.f90 2017-10-23 10:30:22.493326370 +0100 r33406/minpermdist.f90 2017-10-23 10:30:24.061347113 +0100
227: ! NEWMINDIST is always called with PRESERVET .FALSE. here. Hence COORDSA will generally be227: ! NEWMINDIST is always called with PRESERVET .FALSE. here. Hence COORDSA will generally be
228: ! changed to be put into best correspondence with COORDSB.228: ! changed to be put into best correspondence with COORDSB.
229: !229: !
230: ! For TIP potentials specified by W1, W2, W3 and W4 we have to set the atom type correctly230: ! For TIP potentials specified by W1, W2, W3 and W4 we have to set the atom type correctly
231: ! for NEWMINDIST. DJW 23/5/11.231: ! for NEWMINDIST. DJW 23/5/11.
232: !232: !
233:    IF (RBAAT) THEN233:    IF (RBAAT) THEN
234:       CALL RBMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,QBEST,DEBUG)234:       CALL RBMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,QBEST,DEBUG)
235:       CALL QROTMAT (QBEST,RMATBEST)235:       CALL QROTMAT (QBEST,RMATBEST)
236:    ELSEIF (BULKT) THEN236:    ELSEIF (BULKT) THEN
237:       PRINT *,'here A' 
238:          CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
239:       PRINT *,'COORDSA energy=',ENERGY 
240:          CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
241:       PRINT *,'COORDSB energy=',ENERGY 
242:       CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)237:       CALL NEWMINDIST(COORDSB,COORDSA,NATOMS,DISTANCE,BULKT,TWOD,ZSYMSAVE,.FALSE.,RIGIDBODY,DEBUG,RMAT)
243:          CALL POTENTIAL(COORDSA,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
244:       PRINT *,'COORDSA energy=',ENERGY 
245:          CALL POTENTIAL(COORDSB,ENERGY,VNEW,.TRUE.,.FALSE.,RMS,.FALSE.,.FALSE.) 
246:       PRINT *,'COORDSB energy=',ENERGY 
247: !238: !
248: ! How to ensure that we always recognise identical isomers with respect to translation only?239: ! How to ensure that we always recognise identical isomers with respect to translation only?
249: ! Try moving the origin to atom number one.240: ! Try moving the origin to atom number one.
250: !241: !
251: 242: 
252:       XDUMMY=0.0D0243:       XDUMMY=0.0D0
253:       DO J1=1,NATOMS244:       DO J1=1,NATOMS
254:          XDUMMY=XDUMMY + (COORDSA(3*(J1-1)+1)-COORDSB(3*(J1-1)+1)-COORDSA(1)+COORDSB(1) &245:          XDUMMY=XDUMMY + (COORDSA(3*(J1-1)+1)-COORDSB(3*(J1-1)+1)-COORDSA(1)+COORDSB(1) &
255:    &        - BOXLX*NINT((COORDSA(3*(J1-1)+1)-COORDSB(3*(J1-1)+1)-COORDSA(1)+COORDSB(1))/BOXLX))**2 &246:    &        - BOXLX*NINT((COORDSA(3*(J1-1)+1)-COORDSB(3*(J1-1)+1)-COORDSA(1)+COORDSB(1))/BOXLX))**2 &
256:    &                   + (COORDSA(3*(J1-1)+2)-COORDSB(3*(J1-1)+2)-COORDSA(2)+COORDSB(2) &247:    &                   + (COORDSA(3*(J1-1)+2)-COORDSB(3*(J1-1)+2)-COORDSA(2)+COORDSB(2) &


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0