hdiff output

r33431/congrad.f90 2017-11-01 21:30:29.767382466 +0000 r33430/congrad.f90 2017-11-01 21:30:31.107400297 +0000
 10: ! 10: !
 11: !   You should have received a copy of the GNU General Public License 11: !   You should have received a copy of the GNU General Public License
 12: !   along with this program; if not, write to the Free Software 12: !   along with this program; if not, write to the Free Software
 13: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 13: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 14: ! 14: !
 15: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 15: SUBROUTINE CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
 16: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, & 16: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
 17:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, & 17:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
 18:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, JMAXCON, & 18:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,INTIMAGE, KINT, IMSEPMAX, ATOMACTIVE, JMAXCON, &
 19:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, & 19:   &            INTFREEZET, INTFROZEN, CONCUTLOCAL, CONCUT, CONCUTABST, CONCUTABS, CONCUTFRACT, CONCUTFRAC, &
 20:   &  FREEZENODEST, INTSPRINGACTIVET, INTMINFAC, NCONOFF, CONOFFLIST, CONOFFTRIED, INTCONMAX 20:   &  FREEZENODEST, INTSPRINGACTIVET, INTMINFAC
 21: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG 21: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
 22: USE PORFUNCS 22: USE PORFUNCS
 23: IMPLICIT NONE 23: IMPLICIT NONE
 24:             24:            
 25: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX,JMAXNOFF,IMAXNOFF 25: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NREPINT(INTIMAGE+2),ISTAT,NINTMIN,NINTMIN2,MYUNIT,JMAX,IMAX
 26: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX, EMAXNOFF 26: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX
 27: INTEGER JJMAX(INTIMAGE+2) 27: INTEGER JJMAX(INTIMAGE+2)
 28: DOUBLE PRECISION  EEMAX(INTIMAGE+2) 28: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
 29: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1 29: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
 30: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3) 30: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
 31: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI 31: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
 32: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2) 32: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),RMSIM(INTIMAGE+2)
 33: LOGICAL NOINT 33: LOGICAL NOINT
 34: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL 34: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2), CCLOCAL
 35: LOGICAL IMGFREEZE(INTIMAGE) 35: LOGICAL IMGFREEZE(INTIMAGE)
 36: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3) 36: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3)
 44: ECON=0.0D0; EREP=0.0D0 44: ECON=0.0D0; EREP=0.0D0
 45: NINTMIN=0 45: NINTMIN=0
 46: NINTMIN2=0 46: NINTMIN2=0
 47: MYUNIT=6 47: MYUNIT=6
 48:  48: 
 49: EMAX=-1.0D200 49: EMAX=-1.0D200
 50: EEMAX(1:INTIMAGE+2)=-1.0D200 50: EEMAX(1:INTIMAGE+2)=-1.0D200
 51: JJMAX(1:INTIMAGE+2)=-1 51: JJMAX(1:INTIMAGE+2)=-1
 52: JMAX=-1 52: JMAX=-1
 53: IMAX=-1 53: IMAX=-1
 54: EMAXNOFF=-1.0D200 
 55: JMAXNOFF=-1 
 56: IMAXNOFF=-1 
 57:  54: 
 58: ! 55: !
 59: !  Constraint energy and forces. 56: !  Constraint energy and forces.
 60: ! 57: !
 61: DO J2=1,NCONSTRAINT 58: DO J2=1,NCONSTRAINT
 62:    IF (.NOT.CONACTIVE(J2)) CYCLE 59:    IF (.NOT.CONACTIVE(J2)) CYCLE
 63:       CCLOCAL=CONCUTLOCAL(J2) 60:       CCLOCAL=CONCUTLOCAL(J2)
 64:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS 61:       IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
 65:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2) 62:       IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
 66: ! 63: !
 79:       IF (ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL) THEN  76:       IF (ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL) THEN 
 80:          DUMMY=D2-CONDISTREFLOCAL(J2)   77:          DUMMY=D2-CONDISTREFLOCAL(J2)  
 81:          G2(1)=(R2AX-R2BX)/D2 78:          G2(1)=(R2AX-R2BX)/D2
 82:          G2(2)=(R2AY-R2BY)/D2 79:          G2(2)=(R2AY-R2BY)/D2
 83:          G2(3)=(R2AZ-R2BZ)/D2 80:          G2(3)=(R2AZ-R2BZ)/D2
 84:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3) 81:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
 85:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2) 82:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
 86:          EEE(J1)=EEE(J1)  +DUMMY 83:          EEE(J1)=EEE(J1)  +DUMMY
 87:          ECON=ECON        +DUMMY 84:          ECON=ECON        +DUMMY
 88:          IF (DUMMY.GT.EMAX) THEN 85:          IF (DUMMY.GT.EMAX) THEN
 89:             IMAX=J1 86:              IMAX=J1
 90:             JMAX=J2 87:              JMAX=J2
 91:             EMAX=DUMMY 88:              EMAX=DUMMY
 92:          ENDIF 
 93:          IF (DUMMY.GT.EMAXNOFF) THEN 
 94:             IF (.NOT.CONOFFTRIED(J2)) THEN 
 95:                IMAXNOFF=J1 
 96:                JMAXNOFF=J2 
 97:                EMAXNOFF=DUMMY 
 98:             ENDIF 
 99:          ENDIF 89:          ENDIF
100:          IF (DUMMY.GT.EEMAX(J1)) THEN 90:          IF (DUMMY.GT.EEMAX(J1)) THEN
101:             JJMAX(J1)=J2 91:             JJMAX(J1)=J2
102:             EEMAX(J1)=DUMMY 92:             EEMAX(J1)=DUMMY
103:          ENDIF 93:          ENDIF
104:          CONE(J1)=CONE(J1)+DUMMY 94:          CONE(J1)=CONE(J1)+DUMMY
105:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3) 95:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
106:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3) 96:          GGG(NJ1+1:NJ1+3)=GGG(NJ1+1:NJ1+3)-REPGRAD(1:3)
107:       ENDIF 97:       ENDIF
108: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1) 98: !     WRITE(MYUNIT,'(A,2I6,5G20.10)') 'J1,J2,D2,CONDISTREFLOCAL,CCLOCAL,EEE,CONE=',J1,J2,D2,CONDISTREFLOCAL(J2),CCLOCAL,EEE(J1),CONE(J1)
109:    ENDDO 99:    ENDDO
110: ENDDO100: ENDDO
111: IF (JMAX.GT.0) THEN101: IF (JMAX.GT.0) THEN
112:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &102:    WRITE(*,'(A,I6,A,I6,A,2I8)') ' congrad> Highest constraint contribution for any image in image ',IMAX, &
113:  & ' constraint ',JMAX, &103:  & ' constraint ',JMAX, &
114:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF104:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)
115: ENDIF 
116: IF (JMAXNOFF.GT.0) THEN 
117:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, & 
118:  & ' constraint ',JMAXNOFF, & 
119:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF 
120: ELSEIF (JMAX.GT.0) THEN 
121:    JMAXNOFF=JMAX 
122:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad> *** Using highest constraint contribution for any image' 
123:    CONOFFTRIED(1:INTCONMAX)=.FALSE. 
124: ENDIF105: ENDIF
125: JMAXCON=JMAXNOFF106: JMAXCON=JMAX
126: 107: 
127: 108: 
128: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes109: GGG(1:(3*NATOMS))=0.0D0                            ! can delete when loop range above changes
129: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes110: GGG((3*NATOMS)*(INTIMAGE+1)+1:(3*NATOMS)*(INTIMAGE+2))=0.0D0 ! can delete when loop range above changes
130: 111: 
131: ! INTCONST=INTCONSTRAINREPCUT**13112: ! INTCONST=INTCONSTRAINREPCUT**13
132: 113: 
133: EMAX=-1.0D200114: EMAX=-1.0D200
134: EEMAX(1:INTIMAGE+2)=-1.0D200115: EEMAX(1:INTIMAGE+2)=-1.0D200
135: JJMAX(1:INTIMAGE+2)=-1116: JJMAX(1:INTIMAGE+2)=-1
586: 567: 
587: !568: !
588: ! This version of congrad tests for an internal minimum in the569: ! This version of congrad tests for an internal minimum in the
589: ! constraint distances as well as the repulsions.570: ! constraint distances as well as the repulsions.
590: !571: !
591: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)572: SUBROUTINE CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
592: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &573: USE KEY, ONLY: FROZEN, FREEZE, NREPI, NREPJ, NNREPULSIVE, &
593:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &574:   &            NCONSTRAINT, CONI, CONJ, INTCONSTRAINTDEL, CONDISTREF, INTCONSTRAINTREP, CONDISTREFLOCAL, &
594:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &575:   &            CONACTIVE, INTCONSTRAINREPCUT, NREPCUT,FREEZENODEST, INTIMAGE, ATOMACTIVE, KINT, IMSEPMAX, &
595:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &576:   &            INTFREEZET, INTFROZEN, REPI, REPJ, CONCUT, CONCUTLOCAL, &
596:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET, CHECKCONINT, JMAXCON, &577:   &            CONCUTABS, CONCUTABST, CONCUTFRAC, CONCUTFRACT, INTMINFAC, INTSPRINGACTIVET, CHECKCONINT, JMAXCON
597:   &            NCONOFF, CONOFFTRIED, INTCONMAX 
598: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG578: USE COMMONS, ONLY: NATOMS, NOPT, DEBUG
599: IMPLICIT NONE579: IMPLICIT NONE
600:            580:            
601: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2),JMAX,IMAX,JMAXNOFF,IMAXNOFF581: INTEGER :: J1,J2,NI2,NI1,NJ2,NJ1,NMAXINT,NMININT,NCONINT(INTIMAGE+2),NREPINT(INTIMAGE+2),JMAX,IMAX
602: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX, EMAXNOFF582: DOUBLE PRECISION :: ECON, EREP, ETOTAL, RMS, EMAX
603: INTEGER JJMAX(INTIMAGE+2)583: INTEGER JJMAX(INTIMAGE+2)
604: DOUBLE PRECISION  EEMAX(INTIMAGE+2)584: DOUBLE PRECISION  EEMAX(INTIMAGE+2)
605: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1585: DOUBLE PRECISION R1AX,R1AY,R1AZ,R2AX,R2AY,R2AZ,R1BX,R1BY,R1BZ,R2BX,R2BY,R2BZ,D2,D1
606: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)586: DOUBLE PRECISION G1(3),G2(3),DINT,G1INT(3),G2INT(3)
607: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI587: DOUBLE PRECISION DUMMY, REPGRAD(3), INTCONST, D12, DSQ2, DSQ1, DSQI
608: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)588: DOUBLE PRECISION CONE(INTIMAGE+2), REPE(INTIMAGE+2),MAXINT,MININT,REPEINT(INTIMAGE+2),CONEINT(INTIMAGE+2),RMSIMAGE(INTIMAGE+2)
609: LOGICAL NOINT, LPRINT589: LOGICAL NOINT, LPRINT
610: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)590: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), GGG((3*NATOMS)*(INTIMAGE+2)), EEE(INTIMAGE+2)
611: LOGICAL IMGFREEZE(INTIMAGE), PRINTE591: LOGICAL IMGFREEZE(INTIMAGE), PRINTE
612: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), CCLOCAL592: DOUBLE PRECISION DPLUS, ESPRING, SPGRAD(3), CCLOCAL
629: !  Constraint energy and forces.609: !  Constraint energy and forces.
630: !610: !
631: ! For J1 we consider the line segment between image J1-1 and J1.611: ! For J1 we consider the line segment between image J1-1 and J1.
632: ! There are INTIMAGE+1 line segments in total, with an energy contribution612: ! There are INTIMAGE+1 line segments in total, with an energy contribution
633: ! and corresponding gradient terms for each. 613: ! and corresponding gradient terms for each. 
634: ! A and B refer to atoms, 1 and 2 to images J1-1 and J1 corresponding to J1-2 and J1-1 below.614: ! A and B refer to atoms, 1 and 2 to images J1-1 and J1 corresponding to J1-2 and J1-1 below.
635: !615: !
636: ! IMGFREEZE(1:INTIMAGE) refers to the images excluding end points!616: ! IMGFREEZE(1:INTIMAGE) refers to the images excluding end points!
637: !617: !
638: EMAX=-1.0D200618: EMAX=-1.0D200
639: EMAXNOFF=-1.0D200 
640: EEMAX(1:INTIMAGE+2)=-1.0D200619: EEMAX(1:INTIMAGE+2)=-1.0D200
641: JJMAX(1:INTIMAGE+2)=-1620: JJMAX(1:INTIMAGE+2)=-1
642: JMAX=-1621: JMAX=-1
643: IMAX=-1622: IMAX=-1
644: JMAXNOFF=-1 
645: IMAXNOFF=-1 
646: DO J2=1,NCONSTRAINT623: DO J2=1,NCONSTRAINT
647:    IF (.NOT.CONACTIVE(J2)) CYCLE624:    IF (.NOT.CONACTIVE(J2)) CYCLE
648:    CCLOCAL=CONCUTLOCAL(J2)625:    CCLOCAL=CONCUTLOCAL(J2)
649:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS626:    IF (CONCUTABST) CCLOCAL=CCLOCAL+CONCUTABS
650:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)627:    IF (CONCUTFRACT) CCLOCAL=CCLOCAL+CONCUTFRAC*CONDISTREFLOCAL(J2)
651: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG628: !!!!!!!!!!!!!!!!!!!!!!!!!! DEBUG
652:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN629:    IF (INTFROZEN(CONI(J2)).AND.INTFROZEN(CONJ(J2))) THEN
653:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)630:       WRITE(*, '(A,I6,A,2I6)') ' congrad2> ERROR *** constraint ',J2,' between frozen atoms ',CONI(J2),CONJ(J2)
654:       STOP631:       STOP
655:    ENDIF632:    ENDIF
695:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL672:   &       'J1,D2,D1,DINT,MIN diff,CONCUT=',J1,D2,D1,DINT,ABS(D2-CONDISTREFLOCAL(J2)),CCLOCAL
696:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 673:       IF ((ABS(D2-CONDISTREFLOCAL(J2)).GT.CCLOCAL).AND.(J1.LT.INTIMAGE+2)) THEN 
697:          DUMMY=D2-CONDISTREFLOCAL(J2)  674:          DUMMY=D2-CONDISTREFLOCAL(J2)  
698:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)675:          REPGRAD(1:3)=2*INTCONSTRAINTDEL*((DUMMY/CCLOCAL)**2-1.0D0)*DUMMY*G2(1:3)
699:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)676:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
700:          IF (DUMMY.GT.EMAX) THEN677:          IF (DUMMY.GT.EMAX) THEN
701:             IMAX=J1678:             IMAX=J1
702:             JMAX=J2679:             JMAX=J2
703:             EMAX=DUMMY680:             EMAX=DUMMY
704:          ENDIF681:          ENDIF
705:          IF (DUMMY.GT.EMAXNOFF) THEN 
706:             IF (.NOT.CONOFFTRIED(J2)) THEN 
707:                IMAXNOFF=J1 
708:                JMAXNOFF=J2 
709:                EMAXNOFF=DUMMY 
710:             ENDIF 
711:          ENDIF 
712:          IF (DUMMY.GT.EEMAX(J1)) THEN682:          IF (DUMMY.GT.EEMAX(J1)) THEN
713:             JJMAX(J1)=J2683:             JJMAX(J1)=J2
714:             EEMAX(J1)=DUMMY684:             EEMAX(J1)=DUMMY
715:          ENDIF685:          ENDIF
716:          EEE(J1)=EEE(J1)+DUMMY686:          EEE(J1)=EEE(J1)+DUMMY
717:          CONE(J1)=CONE(J1)+DUMMY687:          CONE(J1)=CONE(J1)+DUMMY
718:          ECON=ECON      +DUMMY688:          ECON=ECON      +DUMMY
719:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &689:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'min J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
720:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)690:   &                              SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
721:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)691:          GGG(NI2+1:NI2+3)=GGG(NI2+1:NI2+3)+REPGRAD(1:3)
734:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)704:          DUMMY=INTCONSTRAINTDEL*(DUMMY**2-CCLOCAL**2)**2/(2.0D0*CCLOCAL**2)
735: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN705: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
736: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &706: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'A CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
737: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY707: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
738: !        ENDIF708: !        ENDIF
739:          IF (DUMMY.GT.EMAX) THEN709:          IF (DUMMY.GT.EMAX) THEN
740:             IMAX=J1710:             IMAX=J1
741:             JMAX=J2711:             JMAX=J2
742:             EMAX=DUMMY712:             EMAX=DUMMY
743:          ENDIF713:          ENDIF
744:          IF (DUMMY.GT.EMAXNOFF) THEN 
745:             IF (.NOT.CONOFFTRIED(J2)) THEN 
746:                IMAXNOFF=J1 
747:                JMAXNOFF=J2 
748:                EMAXNOFF=DUMMY 
749:             ENDIF 
750:          ENDIF 
751:          IF (DUMMY.GT.EEMAX(J1-1)) THEN714:          IF (DUMMY.GT.EEMAX(J1-1)) THEN
752:             JJMAX(J1-1)=J2715:             JJMAX(J1-1)=J2
753:             EEMAX(J1-1)=DUMMY716:             EEMAX(J1-1)=DUMMY
754:          ENDIF717:          ENDIF
755:          EEE(J1-1)=EEE(J1-1)+DUMMY718:          EEE(J1-1)=EEE(J1-1)+DUMMY
756:          CONE(J1-1)=CONE(J1-1)+DUMMY719:          CONE(J1-1)=CONE(J1-1)+DUMMY
757:          ECON=ECON      +DUMMY720:          ECON=ECON      +DUMMY
758:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &721:          IF (LPRINT) WRITE(*, '(A,4I6,G15.5)') 'max J1,J2,CONI,CONJ,REPGRAD=',J1,J2,CONI(J2),CONJ(J2), &
759:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)722:   &         SQRT(REPGRAD(1)**2+REPGRAD(2)**2+REPGRAD(3)**2)
760:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)723:          GGG(NI1+1:NI1+3)=GGG(NI1+1:NI1+3)+REPGRAD(1:3)
770: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN733: !        IF (PRINTE.AND.(DUMMY.GT.1.0D-4)) THEN
771: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &734: !           WRITE(*, '(A,2I6,2L5,G20.10)') 'B CONI,CONJ,INTFROZEN(CONI),INTFROZEN(CONJ),DUMMY=', &
772: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY735: ! &                                       CONI(J2),CONJ(J2),INTFROZEN(CONI(J2)),INTFROZEN(CONJ(J2)),DUMMY
773: !        ENDIF736: !        ENDIF
774:          ECON=ECON+DUMMY737:          ECON=ECON+DUMMY
775:          IF (DUMMY.GT.EMAX) THEN738:          IF (DUMMY.GT.EMAX) THEN
776:             IMAX=J1739:             IMAX=J1
777:             JMAX=J2740:             JMAX=J2
778:             EMAX=DUMMY741:             EMAX=DUMMY
779:          ENDIF742:          ENDIF
780:          IF (DUMMY.GT.EMAXNOFF) THEN 
781:             IF (.NOT.CONOFFTRIED(J2)) THEN 
782:                IMAXNOFF=J1 
783:                JMAXNOFF=J2 
784:                EMAXNOFF=DUMMY 
785:             ENDIF 
786:          ENDIF 
787:          IF (DUMMY.GT.EEMAX(J1-1)) THEN743:          IF (DUMMY.GT.EEMAX(J1-1)) THEN
788:             JJMAX(J1-1)=J2744:             JJMAX(J1-1)=J2
789:             EEMAX(J1-1)=DUMMY745:             EEMAX(J1-1)=DUMMY
790:          ENDIF746:          ENDIF
791:          IF (DUMMY.GT.EEMAX(J1)) THEN747:          IF (DUMMY.GT.EEMAX(J1)) THEN
792:             JJMAX(J1)=J2748:             JJMAX(J1)=J2
793:             EEMAX(J1)=DUMMY749:             EEMAX(J1)=DUMMY
794:          ENDIF750:          ENDIF
795:          IF (J1.EQ.2) THEN751:          IF (J1.EQ.2) THEN
796:             EEE(J1)=EEE(J1)+DUMMY752:             EEE(J1)=EEE(J1)+DUMMY
819: !    IF (JJMAX(J2).GT.0) THEN775: !    IF (JJMAX(J2).GT.0) THEN
820: !       WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for image ',J2,' constraint ',JJMAX(J2),' atoms ', &776: !       WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for image ',J2,' constraint ',JJMAX(J2),' atoms ', &
821: !  &                                CONI(JJMAX(J2)),CONJ(JJMAX(J2))777: !  &                                CONI(JJMAX(J2)),CONJ(JJMAX(J2))
822: !    ENDIF778: !    ENDIF
823: ! ENDDO779: ! ENDDO
824: IF (JMAX.GT.0) THEN780: IF (JMAX.GT.0) THEN
825:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for any image in image ',IMAX, &781:    WRITE(*,'(A,I6,A,I6,A,2I6)') ' congrad2> Highest constraint contribution for any image in image ',IMAX, &
826:  & ' constraint ',JMAX, &782:  & ' constraint ',JMAX, &
827:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)783:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX)
828: ENDIF784: ENDIF
829: IF (JMAXNOFF.GT.0) THEN785: JMAXCON=JMAX
830:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> Highest constraint contribution never turned off for any image in image ',IMAXNOFF, & 
831:  & ' constraint ',JMAXNOFF, & 
832:  &                              ' atoms ',CONI(JMAX),CONJ(JMAX),' off=',NCONOFF 
833: ELSEIF (JMAX.GT.0) THEN 
834:    JMAXNOFF=JMAX 
835:    CONOFFTRIED(1:INTCONMAX)=.FALSE. 
836:    WRITE(*,'(A,I6,A,I6,A,2I8,A,I8)') ' congrad2> *** Using highest constraint contribution for any image' 
837: ENDIF 
838: JMAXCON=JMAXNOFF 
839: 786: 
840: ! INTCONST=INTCONSTRAINREPCUT**13787: ! INTCONST=INTCONSTRAINREPCUT**13
841: 788: 
842: EMAX=-1.0D200789: EMAX=-1.0D200
843: EEMAX(1:INTIMAGE+2)=-1.0D200790: EEMAX(1:INTIMAGE+2)=-1.0D200
844: JJMAX(1:INTIMAGE+2)=-1791: JJMAX(1:INTIMAGE+2)=-1
845: JMAX=-1792: JMAX=-1
846: IMAX=-1793: IMAX=-1
847: DO J2=1,NNREPULSIVE794: DO J2=1,NNREPULSIVE
848: !  INTCONST=NREPCUT(J2)**13795: !  INTCONST=NREPCUT(J2)**13


r33431/intlbfgs.f90 2017-11-01 21:30:29.991385447 +0000 r33430/intlbfgs.f90 2017-11-01 21:30:31.331403277 +0000
 12: !   GNU General Public License for more details. 12: !   GNU General Public License for more details.
 13: ! 13: !
 14: !   You should have received a copy of the GNU General Public License 14: !   You should have received a copy of the GNU General Public License
 15: !   along with this program; if not, write to the Free Software 15: !   along with this program; if not, write to the Free Software
 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 16: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 17: ! 17: !
 18: SUBROUTINE INTLBFGS(QSTART,QFINISH) 18: SUBROUTINE INTLBFGS(QSTART,QFINISH)
 19: USE PORFUNCS 19: USE PORFUNCS
 20: USE KEY, ONLY : FREEZENODEST, FREEZETOL, MAXBFGS, CONVR, ATOMSTORES, & 20: USE KEY, ONLY : FREEZENODEST, FREEZETOL, MAXBFGS, CONVR, ATOMSTORES, &
 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, INTMUPDATE, INTDGUESS, & 21:      & INTRMSTOL, INTIMAGE, NREPMAX, NREPULSIVE, INTMUPDATE, INTDGUESS, &
 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, CONOFFLIST, CONOFFTRIED,  & 22:      & NCONSTRAINT, CONI, CONJ, CONDISTREF, INTCONMAX, CONOFFLIST, &
 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, DOBACK, & 33:      & PERMDIST, LOCALPERMCUT, QCILPERMDIST, QCIPDINT, QCIPERMCUT, QCIAMBERT, BONDS, DOBACK, &
 34:      & QCIRESET, QCIRESETINT1, QCIRESETINT2, JMAXCON, NCONOFF 34:      & QCIRESET, QCIRESETINT1, QCIRESETINT2, JMAXCON
 35: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3 35: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3
 36: USE MODCHARMM, ONLY : CHRMMT 36: USE MODCHARMM, ONLY : CHRMMT
 37: USE CHIRALITY 37: USE CHIRALITY
 38:  38: 
 39: IMPLICIT NONE  39: IMPLICIT NONE 
 40:  40: 
 41: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points 41: DOUBLE PRECISION, INTENT(IN) :: QSTART(3*NATOMS), QFINISH(3*NATOMS)  ! The two end points
 42: INTEGER D, U 42: INTEGER D, U
 43: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE, NEIGHBOUR_COORDS(12), CENTRE_COORDS(3) 43: DOUBLE PRECISION DIST, DIST2, RMAT(3,3), SUMEEE, SUMEEE2, SIGMAEEE, NEIGHBOUR_COORDS(12), CENTRE_COORDS(3)
 44: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ 44: DOUBLE PRECISION DMAX, DF, DMIN, LOCALSTEP, ADMAX, DUMMYX, DUMMYY, DUMMYZ
 45: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF 45: INTEGER NDECREASE, NFAIL, NMAXINT, NMININT, JMAX, JMIN, INTIMAGESAVE, NOFF, J1, J2, NQDONE, JA1, JA2, NMOVE, NMOVES, NMOVEF, NCONOFF
 46: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS), STARTGROUP(NPERMGROUP), ENDGROUP(NPERMGROUP) 46: INTEGER PERM(NATOMS), PERMS(NATOMS), PERMF(NATOMS), STARTGROUP(NPERMGROUP), ENDGROUP(NPERMGROUP)
 47: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS), IDENTITY 47: LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM, ADDREP(NATOMS), LDEBUG, REMOVEIMAGE, PERMUTABLE(NATOMS), IDENTITY
 48: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 48: COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 49:  49: 
 50: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY 50: DOUBLE PRECISION DUMMY, DPRAND, DUMMY2, ADUMMY
 51: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF 51: DOUBLE PRECISION BOXLX,BOXLY,BOXLZ,DISTANCE,RMATBEST(3,3),DISTANCES,DISTANCEF
 52: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2 52: INTEGER POINT,NPT,J3,J4,NIMAGEFREEZE,NACTIVE,NBEST,NEWATOM,NBEST2
 53: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND 53: INTEGER TURNONORDER(NATOMS),NBACKTRACK,NQCIFREEZE, NBONDED(NATOMS), BONDEDLIST(NATOMS,6), NBOND
 54: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS), ACID 54: INTEGER NDUMMY, NLASTGOODE, NSTEPSMAX, INGROUP(NATOMS), ACID
 55: LOGICAL CHIRALSR, CHIRALSRP  55: LOGICAL CHIRALSR, CHIRALSRP 
181: ENDIF181: ENDIF
182: IF (INTSTEPS1 < 0) THEN182: IF (INTSTEPS1 < 0) THEN
183:    WRITE(*,'(1x,a)') 'Maximal number of iterations is less than zero! Stop.'183:    WRITE(*,'(1x,a)') 'Maximal number of iterations is less than zero! Stop.'
184:    STOP184:    STOP
185: ENDIF185: ENDIF
186: !186: !
187: ! XYZ, GGG, EEE include the end point images187: ! XYZ, GGG, EEE include the end point images
188: ! X, G do not.188: ! X, G do not.
189: !189: !
190: IF (.NOT.ALLOCATED(CONI)) THEN 190: IF (.NOT.ALLOCATED(CONI)) THEN 
191:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),CONOFFTRIED(INTCONMAX))191:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX))
192:    CONOFFTRIED(1:INTCONMAX)=.FALSE. 
193:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))192:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
194: ENDIF193: ENDIF
195: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))194: X=>XYZ((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
196: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))195: G=>GGG((3*NATOMS)+1:(3*NATOMS)*(INTIMAGE+1))
197: !196: !
198: ! Initialise XYZ197: ! Initialise XYZ
199: !198: !
200: IF (READIMAGET) THEN  ! Note that this will ignore the coordinates in start and finish199: IF (READIMAGET) THEN  ! Note that this will ignore the coordinates in start and finish
201:    LUNIT=GETUNIT()200:    LUNIT=GETUNIT()
202:    OPEN(UNIT=LUNIT,FILE='int.xyz',STATUS='OLD')201:    OPEN(UNIT=LUNIT,FILE='int.xyz',STATUS='OLD')
366:    READ(LUNIT,*) NREPCUT(1:NNREPULSIVE)365:    READ(LUNIT,*) NREPCUT(1:NNREPULSIVE)
367:    WRITE(*,'(A)') ' intlbfgs> read NREPCUT:'366:    WRITE(*,'(A)') ' intlbfgs> read NREPCUT:'
368: !  WRITE(*,'(6G20.10)') NREPCUT(1:NNREPULSIVE)367: !  WRITE(*,'(6G20.10)') NREPCUT(1:NNREPULSIVE)
369:    READ(LUNIT,*) INTFROZEN(1:NATOMS)368:    READ(LUNIT,*) INTFROZEN(1:NATOMS)
370:    WRITE(*,'(A)') ' intlbfgs> read INTFROZEN'369:    WRITE(*,'(A)') ' intlbfgs> read INTFROZEN'
371: !  WRITE(*,'(12L5)') INTFROZEN(1:NATOMS)370: !  WRITE(*,'(12L5)') INTFROZEN(1:NATOMS)
372: 371: 
373:    NCONOFF=0372:    NCONOFF=0
374:    READ(LUNIT,*,END=742) NCONOFF373:    READ(LUNIT,*,END=742) NCONOFF
375:    IF (NCONOFF.GT.0) READ(LUNIT,*) CONOFFLIST(1:NCONOFF)374:    IF (NCONOFF.GT.0) READ(LUNIT,*) CONOFFLIST(1:NCONOFF)
376:    IF (NCONOFF.GT.0) READ(LUNIT,*) CONOFFTRIED(1:NCONOFF) 
377: 742 CONTINUE375: 742 CONTINUE
378:    CLOSE(LUNIT)376:    CLOSE(LUNIT)
379: 377: 
380:    GLAST(1:D)=G(1:D)378:    GLAST(1:D)=G(1:D)
381:    XSAVE(1:D)=X(1:D)379:    XSAVE(1:D)=X(1:D)
382:    GOTO 986380:    GOTO 986
383: ENDIF381: ENDIF
384: IF (INTFREEZET) THEN382: IF (INTFREEZET) THEN
385:    DO J1=1,NATOMS383:    DO J1=1,NATOMS
386:       IF (INTFROZEN(J1)) THEN384:       IF (INTFROZEN(J1)) THEN
587: ! ELSE585: ! ELSE
588:    ! CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)586:    ! CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
589: ENDIF587: ENDIF
590: EOLD=ETOTAL588: EOLD=ETOTAL
591: GLAST(1:D)=G(1:D)589: GLAST(1:D)=G(1:D)
592: XSAVE(1:D)=X(1:D)590: XSAVE(1:D)=X(1:D)
593: 591: 
594: IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN592: IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN
595:    WRITE(*,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=', &593:    WRITE(*,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=', &
596:   &                       ETOTAL/INTIMAGE,COLDFUSIONLIMIT594:   &                       ETOTAL/INTIMAGE,COLDFUSIONLIMIT
597:    DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST,CONOFFTRIED)595:    DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST)
598:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &596:    DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
599:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)597:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
600:    INTIMAGE=INTIMAGESAVE598:    INTIMAGE=INTIMAGESAVE
601:    RETURN599:    RETURN
602: ENDIF600: ENDIF
603: 601: 
604: ! IF (DEBUG) WRITE(*,'(A6,A20,A20,A9,A9)') 'Iter','Energy per image','RMS Force','Step'602: ! IF (DEBUG) WRITE(*,'(A6,A20,A20,A9,A9)') 'Iter','Energy per image','RMS Force','Step'
605: 603: 
606: !604: !
607: ! In this block PERMGROUP(NDUMMY+J2-1) counts through the atom indices specified in perm.allow,605: ! In this block PERMGROUP(NDUMMY+J2-1) counts through the atom indices specified in perm.allow,
625: 623: 
626: 567 CONTINUE624: 567 CONTINUE
627: 625: 
628: DO ! Main do loop with counter NITERDONE, initially set to one626: DO ! Main do loop with counter NITERDONE, initially set to one
629: 627: 
630: !628: !
631: ! Are we stuck? If so, try resetting problem atoms to previous image.629: ! Are we stuck? If so, try resetting problem atoms to previous image.
632: !630: !
633: IF (QCIRESET) THEN631: IF (QCIRESET) THEN
634: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN632: !  IF ((SWITCHED.AND.(MOD(NITERDONE-1,QCIRESETINT2).EQ.0)).OR.((.NOT.SWITCHED).AND.(MOD(NITERDONE-1,QCIRESETINT1).EQ.0))) THEN
635: !  PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1633:    PRINT *,'intlbfgs> NITERDONE,NLASTGOODE,QCIRESETINT1=',NITERDONE,NLASTGOODE,QCIRESETINT1
636:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN634:    IF ((.NOT.SWITCHED).AND.(NITERDONE-NLASTGOODE.GT.QCIRESETINT1)) THEN
637:       CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)635:       CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
638:       CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)636:       CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
639:       WRITE(*,'(2(A,I6))') 'intlbfgs> Interpolation seems to be stuck. Turn off worst constraint ',JMAXCON,' total off=',NCONOFF+1637:       WRITE(*,'(A,I6)') 'intlbfgs> Interpolation seems to be stuck. Turn off worst constraint ',JMAXCON
640:       IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN638:       IF ((JMAXCON.LT.1).OR.(JMAXCON.GT.NCONSTRAINT)) THEN
641:          WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'639:          WRITE(*,'(A)') 'intlbfgs> *** ERROR *** constraint index out of allowed range'
642:          WRITE(*,'(A,I6)') 'NCONSTRAINT,NCONOFF=',NCONOFF 
643:          WRITE(*,'(A)') 'CONOFFTRIED:' 
644:          WRITE(*,'(20L5)') CONOFFTRIED(1:NCONSTRAINT) 
645:          STOP640:          STOP
646:       ENDIF641:       ENDIF
647:       NCONOFF=NCONOFF+1642:       NCONOFF=NCONOFF+1
648:       CONOFFLIST(NCONOFF)=JMAXCON643:       CONOFFLIST(NCONOFF)=JMAXCON
649:       CONACTIVE(JMAXCON)=.FALSE.644:       CONACTIVE(JMAXCON)=.FALSE.
650:       CONOFFTRIED(JMAXCON)=.TRUE. 
651:       NLASTGOODE=NITERDONE645:       NLASTGOODE=NITERDONE
652:       LASTGOODE=ETOTAL646:       LASTGOODE=ETOTAL
653: !     STOP647: !     STOP
654:    ENDIF648:    ENDIF
655: ENDIF649: ENDIF
656: 650: 
657: !651: !
658: !  Check permutational alignments. Maintain a list of the permutable groups where all652: !  Check permutational alignments. Maintain a list of the permutable groups where all
659: !  members are active. See if we have any new complete groups. MUST update NDUMMY653: !  members are active. See if we have any new complete groups. MUST update NDUMMY
660: !  counter to step through permutable atom list.654: !  counter to step through permutable atom list.
916: !!!910: !!!
917: !!! Use the minimum of the end point distances and INTCONSTRAINREPCUT for each contact.911: !!! Use the minimum of the end point distances and INTCONSTRAINREPCUT for each contact.
918: !!!912: !!!
919: !!                     REPCUT(NREPULSIVE)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)913: !!                     REPCUT(NREPULSIVE)=MIN(DMIN-1.0D-3,INTCONSTRAINREPCUT)
920: !!                  ENDDO myrep2914: !!                  ENDDO myrep2
921: !!               ENDDO915: !!               ENDDO
922: !!               WRITE(*,'(A,I6,A)') ' intlbfgs> Now it looks like there are ',NREPULSIVE,' possible repulsions before adding new atom'916: !!               WRITE(*,'(A,I6,A)') ' intlbfgs> Now it looks like there are ',NREPULSIVE,' possible repulsions before adding new atom'
923: !!!!!!!!!!!!!!!DEBUG DJW !!!!!!!!!!!917: !!!!!!!!!!!!!!!DEBUG DJW !!!!!!!!!!!
924: 918: 
925:       IF (NCONOFF.GT.0) THEN919:       IF (NCONOFF.GT.0) THEN
926:          CONACTIVE(CONOFFLIST(NCONOFF))=.TRUE.920:          WRITE(*,'(A,I6)') 'intlbfgs> Turn back on constraint ',CONOFFLIST(NCONOFF)
927:          WRITE(*,'(2(A,I6))') 'intlbfgs> Turn back on constraint ',CONOFFLIST(NCONOFF),' total off=',NCONOFF-1921:          CONACTIVE(NCONOFF)=.TRUE.
928:          NCONOFF=NCONOFF-1922:          NCONOFF=NCONOFF-1
929:       ELSE923:       ELSE
930:          CALL DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE,AABACK,BACKDONE)  924:          CALL DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE,AABACK,BACKDONE)  
931:       ENDIF925:       ENDIF
932:       NLASTGOODE=NITERDONE926:       NLASTGOODE=NITERDONE
933:       LASTGOODE=ETOTAL927:       LASTGOODE=ETOTAL
934:    ENDIF928:    ENDIF
935:    GTMP(1:D)=0.0D0929:    GTMP(1:D)=0.0D0
936:    CALL MAKESTEP(NITERUSE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)930:    CALL MAKESTEP(NITERUSE,POINT,DIAG,INTIMAGE,SEARCHSTEP,G,GTMP,STP,GDIF,NPT,D,RHO1,ALPHA)
937: !931: !
1338:       RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))1332:       RMS=SQRT(RMS/((3*NATOMS)*INTIMAGE))
1339:       EEE(1:INTIMAGE+2)=USEFRAC*EEETMP(1:INTIMAGE+2)+(1.0D0-USEFRAC)*EEE(1:INTIMAGE+2)1333:       EEE(1:INTIMAGE+2)=USEFRAC*EEETMP(1:INTIMAGE+2)+(1.0D0-USEFRAC)*EEE(1:INTIMAGE+2)
1340:       WORST=-1.0D1001334:       WORST=-1.0D100
1341:       DO J4=2,INTIMAGE+11335:       DO J4=2,INTIMAGE+1
1342:          IF (EEE(J4).GT.WORST) WORST=EEE(J4)1336:          IF (EEE(J4).GT.WORST) WORST=EEE(J4)
1343:       ENDDO1337:       ENDDO
1344:       IF (DEBUG) WRITE(*,'(A,G20.10,A,I8)') 'intlbfgs> Highest QCI image energy=',WORST,' images=',INTIMAGE1338:       IF (DEBUG) WRITE(*,'(A,G20.10,A,I8)') 'intlbfgs> Highest QCI image energy=',WORST,' images=',INTIMAGE
1345:    ENDIF1339:    ENDIF
1346:    IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN1340:    IF (ETOTAL/INTIMAGE.LT.COLDFUSIONLIMIT) THEN
1347:       WRITE(*,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/INTIMAGE,COLDFUSIONLIMIT1341:       WRITE(*,'(A,2G20.10)') ' intlbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/INTIMAGE,COLDFUSIONLIMIT
1348:       DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST,CONOFFTRIED)1342:       DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST)
1349:       DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &1343:       DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
1350:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)1344:   &              DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
1351:       QCIIMAGE=INTIMAGE1345:       QCIIMAGE=INTIMAGE
1352:       INTIMAGE=INTIMAGESAVE1346:       INTIMAGE=INTIMAGESAVE
1353:       RETURN1347:       RETURN
1354:    ENDIF1348:    ENDIF
1355: 1349: 
1356:    STEPTOT = SUM(STEPIMAGE)/INTIMAGE1350:    STEPTOT = SUM(STEPIMAGE)/INTIMAGE
1357: 1351: 
1358:    MAXRMS=-1.0D01352:    MAXRMS=-1.0D0
1512:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW1506:    INTDGUESS=DIAG(1) ! should be ok for subsequent runs of the same system DJW
1513:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 1507:    IF ((.NOT.SWITCHED).AND.(MAXRMS<=INTRMSTOL).AND.NITERDONE>1) EXITSTATUS=1 
1514:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 1508:    IF (SWITCHED.AND.(MAXRMS<=CONVR).AND.NITERDONE>1) EXITSTATUS=1 
1515:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=21509:    IF (NITERDONE==NSTEPSMAX) EXITSTATUS=2
1516:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW1510:    IF ((.NOT.SWITCHED).AND.(MOD(NITERDONE,INTRELSTEPS).EQ.0)) EXITSTATUS=1 ! Add an atom every INTRELSTEPS !!! DJW
1517: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX1511: !  PRINT '(A,2G20.10,3I8)','MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX=',MAXRMS,INTRMSTOL,NITERDONE,NITERDONE,NSTEPSMAX
1518: 1512: 
1519:    IF (EXITSTATUS > 0) THEN  1513:    IF (EXITSTATUS > 0) THEN  
1520:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on1514:       IF ((.NOT.SWITCHED).AND.(EXITSTATUS.EQ.1)) THEN ! add active atom or restart with true potential on
1521: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771515: !        IF (ETOTAL/INTIMAGE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1522:          PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))1516:          PRINT '(A,3G20.10)','MAXEEE,MAXCONE,scaled=',MAXEEE,MAXCONE,MAXCONE*MAX(0.2D0,NACTIVE*1.0D0/(NATOMS*1.0D0))
1523:          IF (MAXEEE.GT.MAXCONE*MAX(0.1D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 7771517:          IF (MAXEEE.GT.MAXCONE*MAX(0.2D0,NACTIVE*1.0D0/(NATOMS*1.0D0))) GOTO 777
1524:          IF (NACTIVE.LT.NATOMS) THEN 1518:          IF (NACTIVE.LT.NATOMS) THEN 
1525:             ADDATOM=.TRUE.1519:             ADDATOM=.TRUE.
1526:             GOTO 7771520:             GOTO 777
1527:          ENDIF1521:          ENDIF
1528:          CALL MYCPU_TIME(FTIME,.FALSE.)1522:          CALL MYCPU_TIME(FTIME,.FALSE.)
1529:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &1523:          WRITE(*,'(A,I6,A,F12.6,A,I6,A,G20.10)') ' intlbfgs> switch on true potential at step ',NITERDONE, &
1530:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME1524:   &                                     ' fraction=',INTCONFRAC,' images=',INTIMAGE,' time=',FTIME-STIME
1531:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1525:          IF (DEBUG) CALL INTRWG(NACTIVE,NITERDONE,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1532:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)1526:          IF (DEBUG) CALL WRITEPROFILE(NITERDONE,EEE,INTIMAGE)
1533:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'1527:          WRITE(*,'(A,I6,A,F15.6)') ' intlbfgs> Allowing ',INTCONSTEPS,' further optimization steps'
1615: BESTWORST=WORST1609: BESTWORST=WORST
1616: BESTINTIMAGE=INTIMAGE1610: BESTINTIMAGE=INTIMAGE
1617: IF (ALLOCATED(QCIXYZ)) DEALLOCATE(QCIXYZ)1611: IF (ALLOCATED(QCIXYZ)) DEALLOCATE(QCIXYZ)
1618: ALLOCATE(QCIXYZ(3*NATOMS*(INTIMAGE+2)))1612: ALLOCATE(QCIXYZ(3*NATOMS*(INTIMAGE+2)))
1619: QCIXYZ(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))1613: QCIXYZ(1:3*NATOMS*(INTIMAGE+2))=XYZ(1:3*NATOMS*(INTIMAGE+2))
1620: WRITE(*,'(A,I8,A,G20.10)') 'intlbfgs> retaining ',INTIMAGE,' QCI images, highest energy=',BESTWORST1614: WRITE(*,'(A,I8,A,G20.10)') 'intlbfgs> retaining ',INTIMAGE,' QCI images, highest energy=',BESTWORST
1621: 1615: 
1622: CALL INTRWG(NACTIVE,0,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1616: CALL INTRWG(NACTIVE,0,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1623: CALL WRITEPROFILE(0,EEE,INTIMAGE)1617: CALL WRITEPROFILE(0,EEE,INTIMAGE)
1624: 1618: 
1625: DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST,CONOFFTRIED)1619: DEALLOCATE(CONI,CONJ,CONDISTREF,REPI,REPJ,NREPI,NREPJ,REPCUT,NREPCUT,CONCUT,CONOFFLIST)
1626: DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &1620: DEALLOCATE(TRUEEE, EEETMP, MYGTMP, GTMP, &
1627:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)1621:   &      DIAG, STP, SEARCHSTEP, GDIF,GLAST, XSAVE, XYZ, GGG, CHECKG, IMGFREEZE, EEE, STEPIMAGE)
1628: QCIIMAGE=INTIMAGE1622: QCIIMAGE=INTIMAGE
1629: INTIMAGE=INTIMAGESAVE1623: INTIMAGE=INTIMAGESAVE
1630: 1624: 
1631: END SUBROUTINE INTLBFGS1625: END SUBROUTINE INTLBFGS
1632: !1626: !
1633: ! Neighbour list for repulsions to reduce cost of constraint potential.1627: ! Neighbour list for repulsions to reduce cost of constraint potential.
1634: !1628: !
1635: SUBROUTINE CHECKREP(INTIMAGE,XYZ,NOPT,NNSTART,NSTART)1629: SUBROUTINE CHECKREP(INTIMAGE,XYZ,NOPT,NNSTART,NSTART)
1696:    ENDDO 1690:    ENDDO 
1697: 246 CONTINUE1691: 246 CONTINUE
1698: ENDDO1692: ENDDO
1699: IF (DEBUG) WRITE(*,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE1693: IF (DEBUG) WRITE(*,'(A,2I8)') ' checkrep> number of active repulsions and total=',NNREPULSIVE,NREPULSIVE
1700: 1694: 
1701: END SUBROUTINE CHECKREP1695: END SUBROUTINE CHECKREP
1702: 1696: 
1703: SUBROUTINE INTRWG(NACTIVE,NITER,INTIMAGE,XYZ,TURNONORDER,NCONOFF)1697: SUBROUTINE INTRWG(NACTIVE,NITER,INTIMAGE,XYZ,TURNONORDER,NCONOFF)
1704: USE PORFUNCS1698: USE PORFUNCS
1705: USE KEY,ONLY: STOCKT,STOCKAAT, RBAAT, ATOMACTIVE, NCONSTRAINT, CONACTIVE, NREPULSIVE, NNREPULSIVE, REPI, REPJ, REPCUT, NREPCUT, &1699: USE KEY,ONLY: STOCKT,STOCKAAT, RBAAT, ATOMACTIVE, NCONSTRAINT, CONACTIVE, NREPULSIVE, NNREPULSIVE, REPI, REPJ, REPCUT, NREPCUT, &
1706:   &           NREPMAX, NREPI, NREPJ, INTFROZEN, CONOFFLIST,CONOFFTRIED1700:   &           NREPMAX, NREPI, NREPJ, INTFROZEN, CONOFFLIST
1707: USE COMMONS, ONLY: NATOMS, DEBUG1701: USE COMMONS, ONLY: NATOMS
1708: IMPLICIT NONE1702: IMPLICIT NONE
1709: INTEGER NCONOFF1703: INTEGER NCONOFF
1710: CHARACTER(LEN=10) :: XYZFILE   = 'int.xyz   '1704: CHARACTER(LEN=10) :: XYZFILE   = 'int.xyz   '
1711: CHARACTER(LEN=10) :: QCIFILE   = 'QCIdump   '1705: CHARACTER(LEN=10) :: QCIFILE   = 'QCIdump   '
1712: INTEGER,INTENT(IN) :: NITER, TURNONORDER(NATOMS)1706: INTEGER,INTENT(IN) :: NITER, TURNONORDER(NATOMS)
1713: INTEGER :: J1,J2,INTIMAGE,J3,NACTIVE,LUNIT,GETUNIT1707: INTEGER :: J1,J2,INTIMAGE,J3,NACTIVE,LUNIT,GETUNIT
1714: CHARACTER(LEN=80) :: FILENAME,DUMMYS1708: CHARACTER(LEN=80) :: FILENAME,DUMMYS
1715: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2))1709: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2))
1716: 1710: 
1717: FILENAME=XYZFILE1711: FILENAME=XYZFILE
1737: ENDDO1731: ENDDO
1738: 1732: 
1739: WRITE(*,*) 'rwg> Interpolated image coordinates were saved to xyz file "'//TRIM(FILENAME)//'"'1733: WRITE(*,*) 'rwg> Interpolated image coordinates were saved to xyz file "'//TRIM(FILENAME)//'"'
1740: 1734: 
1741: CLOSE(LUNIT)1735: CLOSE(LUNIT)
1742: 1736: 
1743: FILENAME=QCIFILE1737: FILENAME=QCIFILE
1744: LUNIT=GETUNIT()1738: LUNIT=GETUNIT()
1745: OPEN(UNIT=LUNIT,FILE=TRIM(ADJUSTL(FILENAME)),STATUS='replace')1739: OPEN(UNIT=LUNIT,FILE=TRIM(ADJUSTL(FILENAME)),STATUS='replace')
1746: 1740: 
1747: IF (DEBUG) WRITE(*,'(A,I10,A)') ' intlbfgs> dumping state for ',NACTIVE,' active atoms'1741: WRITE(*,'(A,I10,A)') ' intlbfgs> dumping state for ',NACTIVE,' active atoms'
1748: WRITE(LUNIT,'(I10)') NACTIVE1742: WRITE(LUNIT,'(I10)') NACTIVE
1749: ! WRITE(*,'(A,I10,A)') ' intlbfgs> dumping turnonorder for ',NACTIVE,' active atoms'1743: WRITE(*,'(A,I10,A)') ' intlbfgs> dumping turnonorder for ',NACTIVE,' active atoms'
1750: WRITE(LUNIT,'(12I8)') TURNONORDER(1:NACTIVE)1744: WRITE(LUNIT,'(12I8)') TURNONORDER(1:NACTIVE)
1751: ! WRITE(*,'(A)') ' intlbfgs> dumping atomactive'1745: WRITE(*,'(A)') ' intlbfgs> dumping atomactive'
1752: WRITE(LUNIT,'(12L5)') ATOMACTIVE(1:NATOMS)1746: WRITE(LUNIT,'(12L5)') ATOMACTIVE(1:NATOMS)
1753: WRITE(LUNIT,'(I10)') NCONSTRAINT1747: WRITE(LUNIT,'(I10)') NCONSTRAINT
1754: ! WRITE(*,'(A,I10,A)') ' intlbfgs> dumping conactive for ',NCONSTRAINT,' constraints'1748: WRITE(*,'(A,I10,A)') ' intlbfgs> dumping conactive for ',NCONSTRAINT,' constraints'
1755: WRITE(LUNIT,'(12L5)') CONACTIVE(1:NCONSTRAINT)1749: WRITE(LUNIT,'(12L5)') CONACTIVE(1:NCONSTRAINT)
1756: 1750: 
1757:    WRITE(LUNIT,'(3I12,G20.10)') NREPULSIVE,NNREPULSIVE,NREPMAX1751:    WRITE(LUNIT,'(3I12,G20.10)') NREPULSIVE,NNREPULSIVE,NREPMAX
1758:    ! WRITE(*,'(A,3I10,G20.10)') 'intlbfgs> dumping NREPULSIVE,NNREPULSIVE,NREPMAX=',NREPULSIVE,NNREPULSIVE,NREPMAX1752:    WRITE(*,'(A,3I10,G20.10)') 'intlbfgs> dumping NREPULSIVE,NNREPULSIVE,NREPMAX=',NREPULSIVE,NNREPULSIVE,NREPMAX
1759: 1753: 
1760:    WRITE(LUNIT,'(12I8)') REPI(1:NREPULSIVE)1754:    WRITE(LUNIT,'(12I8)') REPI(1:NREPULSIVE)
1761:    ! WRITE(*,'(A)') ' intlbfgs> dumped REPI:'1755:    WRITE(*,'(A)') ' intlbfgs> dumped REPI:'
1762:    WRITE(LUNIT,'(12I8)') REPJ(1:NREPULSIVE)1756:    WRITE(LUNIT,'(12I8)') REPJ(1:NREPULSIVE)
1763:    ! WRITE(*,'(A)') ' intlbfgs> dumped REPJ:'1757:    WRITE(*,'(A)') ' intlbfgs> dumped REPJ:'
1764:    WRITE(LUNIT,'(12I8)') NREPI(1:NNREPULSIVE)1758:    WRITE(LUNIT,'(12I8)') NREPI(1:NNREPULSIVE)
1765:    ! WRITE(*,'(A)') ' intlbfgs> dumped NREPI:'1759:    WRITE(*,'(A)') ' intlbfgs> dumped NREPI:'
1766:    WRITE(LUNIT,'(12I8)') NREPJ(1:NNREPULSIVE)1760:    WRITE(LUNIT,'(12I8)') NREPJ(1:NNREPULSIVE)
1767:    ! WRITE(*,'(A)') ' intlbfgs> dumped NREPJ:'1761:    WRITE(*,'(A)') ' intlbfgs> dumped NREPJ:'
1768: 1762: 
1769:    WRITE(LUNIT,'(6G20.10)') REPCUT(1:NREPULSIVE)1763:    WRITE(LUNIT,'(6G20.10)') REPCUT(1:NREPULSIVE)
1770:    ! WRITE(*,'(A)') ' intlbfgs> dumped REPCUT:'1764:    WRITE(*,'(A)') ' intlbfgs> dumped REPCUT:'
1771:    WRITE(LUNIT,'(6G20.10)') NREPCUT(1:NNREPULSIVE)1765:    WRITE(LUNIT,'(6G20.10)') NREPCUT(1:NNREPULSIVE)
1772:    ! WRITE(*,'(A)') ' intlbfgs> dumped NREPCUT:'1766:    WRITE(*,'(A)') ' intlbfgs> dumped NREPCUT:'
1773: 1767: 
1774:    WRITE(LUNIT,'(12L5)') INTFROZEN(1:NATOMS)1768:    WRITE(LUNIT,'(12L5)') INTFROZEN(1:NATOMS)
1775:    ! WRITE(*,'(A)') ' intlbfgs> dumped INTFROZEN'1769:    WRITE(*,'(A)') ' intlbfgs> dumped INTFROZEN'
1776: 1770: 
1777:    WRITE(LUNIT,'(I8)') NCONOFF1771:    WRITE(LUNIT,'(I8)') NCONOFF
1778:    IF (NCONOFF.GT.0) WRITE(LUNIT,'(12I8)') CONOFFLIST(1:NCONOFF)1772:    IF (NCONOFF.GT.0) WRITE(LUNIT,'(12I8)') CONOFFLIST(1:NCONOFF)
1779:    IF (NCONOFF.GT.0) WRITE(LUNIT,'(12L5)') CONOFFTRIED(1:NCONOFF)1773:    WRITE(*,'(A)') ' intlbfgs> dumped NCONOFF and CONOFFLIST'
1780:    ! WRITE(*,'(A)') ' intlbfgs> dumped NCONOFF and CONOFFLIST' 
1781: 1774: 
1782: CLOSE(LUNIT)1775: CLOSE(LUNIT)
1783: 1776: 
1784: END SUBROUTINE INTRWG1777: END SUBROUTINE INTRWG
1785: 1778: 
1786: SUBROUTINE WRITEPROFILE(NITER,EEE,INTIMAGE)1779: SUBROUTINE WRITEPROFILE(NITER,EEE,INTIMAGE)
1787: IMPLICIT NONE 1780: IMPLICIT NONE 
1788: INTEGER,INTENT(IN) :: NITER, INTIMAGE1781: INTEGER,INTENT(IN) :: NITER, INTIMAGE
1789: INTEGER :: I,LUNIT,GETUNIT1782: INTEGER :: I,LUNIT,GETUNIT
1790: DOUBLE PRECISION :: EEE(INTIMAGE+2)1783: DOUBLE PRECISION :: EEE(INTIMAGE+2)
1808: CLOSE(LUNIT)1801: CLOSE(LUNIT)
1809: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'1802: WRITE(*,'(A)') ' writeprofile> Interpolated energy profile was saved to file "'//trim(filename)//'"'
1810: 1803: 
1811: END SUBROUTINE WRITEPROFILE1804: END SUBROUTINE WRITEPROFILE
1812: 1805: 
1813: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE,AABACK,BACKDONE)1806: SUBROUTINE DOADDATOM(NCONSTRAINT,NTRIES,NEWATOM,IMGFREEZE,INTIMAGE,XYZ,EEE,GGG,TURNONORDER,NITERDONE,NACTIVE,AABACK,BACKDONE)
1814: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &1807: USE KEY, ONLY : CONACTIVE, CONI, CONJ, ATOMACTIVE, CONDISTREF, REPI, REPJ, REPCUT, INTREPSEP,  &
1815:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &1808:   &             INTCONSTRAINREPCUT, NREPULSIVE, NREPMAX, MAXCONUSE, CHECKCONINT, &
1816:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, ATOMSTORES, QCIADDACIDT, &1809:   &             FREEZENODEST, NNREPULSIVE, INTFROZEN, ATOMSTORES, QCIADDACIDT, &
1817:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &1810:   &             NREPULSIVEFIX, REPIFIX, REPJFIX, REPCUTFIX, NREPI, NREPJ, NREPCUT, MAXNACTIVE, &
1818:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, &1811:   &             NCONSTRAINTFIX, CONIFIX, CONJFIX, INTCONCUT, INTCONSEP, QCIRADSHIFTT, QCIRADSHIFT, QCIADDREP, DOBACK, DOBACKALL
1819:   &             CONOFFTRIED, QCIADDREP, DOBACK, DOBACKALL 
1820: USE COMMONS, ONLY: NATOMS, DEBUG1812: USE COMMONS, ONLY: NATOMS, DEBUG
1821: IMPLICIT NONE1813: IMPLICIT NONE
1822: INTEGER INTIMAGE1814: INTEGER INTIMAGE
1823: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &1815: INTEGER NBEST, NCONTOACTIVE(NATOMS),  NCONSTRAINT, J2, NTRIES(NATOMS), NEWATOM,  CONLIST(NATOMS), N1, N2, N3, &
1824:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE1816:   &     NTOADD, NADDED, NMININT, NMAXINT, TURNONORDER(NATOMS), NDUMMY, J1, J3, NITERDONE, NCONFORNEWATOM, NACTIVE
1825: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN1817: DOUBLE PRECISION DUMMY, DUMMY2, DPRAND, RANDOM, CONDIST(NATOMS), DMIN
1826: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID1818: INTEGER NDFORNEWATOM, BESTPRESERVEDN(NATOMS), ACID
1827: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)1819: DOUBLE PRECISION BESTPRESERVEDD(NATOMS), BESTCLOSESTD(NATOMS), INVDTOACTIVE(NATOMS)
1828: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID, AABACK(NATOMS), BACKDONE, FTEST1820: LOGICAL IMGFREEZE(INTIMAGE), ADDREP(NATOMS), CHOSENACID, AABACK(NATOMS), BACKDONE
1829: DOUBLE PRECISION P1(3), P2(3), P3(3), SOL1(3), SOL2(3), R1, R2, R3 
1830: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE01821: DOUBLE PRECISION C1, C2, C3, VEC1(3), VEC2(3), VEC3(3), ESAVED, ESAVEC, ESAVE0
1831: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE1822: INTEGER NCFORNEWATOM, BESTCLOSESTN(NATOMS), NNREPSAVE, NREPSAVE
1832: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), XSAVED(3,INTIMAGE+2), XSAVEC(3,INTIMAGE+2), XSAVE0(3,INTIMAGE+2),FRAC,RAN1, &1823: DOUBLE PRECISION XYZ((3*NATOMS)*(INTIMAGE+2)), XSAVED(3,INTIMAGE+2), XSAVEC(3,INTIMAGE+2), XSAVE0(3,INTIMAGE+2),FRAC,RAN1, &
1833:   &              RMS,EEE(INTIMAGE+2),GGG((3*NATOMS)*(INTIMAGE+2)),ETOTAL,DS,DF,DNORM,D1SQ,D2SQ1824:   &              RMS,EEE(INTIMAGE+2),GGG((3*NATOMS)*(INTIMAGE+2)),ETOTAL,DS,DF,DNORM
1834: 1825: 
1835: 1826: 
1836: NTOADD=11827: NTOADD=1
1837: NADDED=01828: NADDED=0
1838: CHOSENACID=.FALSE.1829: CHOSENACID=.FALSE.
1839: 1830: 
1840: IF (DOBACK.AND.(.NOT.BACKDONE)) THEN1831: IF (DOBACK.AND.(.NOT.BACKDONE)) THEN
1841:    DO J1=1,NATOMS1832:    DO J1=1,NATOMS
1842:       IF (AABACK(J1)) THEN1833:       IF (AABACK(J1)) THEN
1843:          IF (.NOT.ATOMACTIVE(J1)) GOTO 7631834:          IF (.NOT.ATOMACTIVE(J1)) GOTO 763
1861: ! calls to CHECKREP1852: ! calls to CHECKREP
1862: !1853: !
1863: NNREPSAVE=NNREPULSIVE1854: NNREPSAVE=NNREPULSIVE
1864: NREPSAVE=NREPULSIVE1855: NREPSAVE=NREPULSIVE
1865: 542   CONTINUE1856: 542   CONTINUE
1866: !     DUMMY=1.0D1001857: !     DUMMY=1.0D100
1867:       NBEST=01858:       NBEST=0
1868:       NCONTOACTIVE(1:NATOMS)=01859:       NCONTOACTIVE(1:NATOMS)=0
1869:       INVDTOACTIVE(1:NATOMS)=0.0D01860:       INVDTOACTIVE(1:NATOMS)=0.0D0
1870:       DO J2=1,NCONSTRAINT1861:       DO J2=1,NCONSTRAINT
1871:          IF (CONACTIVE(J2)) CYCLE     ! count new, inactive constraints1862:          IF (CONACTIVE(J2)) CYCLE   ! count new, inactive constraints
1872:          IF (CONOFFTRIED(J2)) CYCLE   ! if we've tried turning it off, it must actually be active. Don't try again. 
1873:          IF (ATOMACTIVE(CONI(J2))) THEN1863:          IF (ATOMACTIVE(CONI(J2))) THEN
1874:             IF (.NOT.ATOMACTIVE(CONJ(J2))) THEN1864:             IF (.NOT.ATOMACTIVE(CONJ(J2))) THEN
1875:                NCONTOACTIVE(CONJ(J2))=NCONTOACTIVE(CONJ(J2))+11865:                NCONTOACTIVE(CONJ(J2))=NCONTOACTIVE(CONJ(J2))+1
1876:                IF (1.0D0/CONDISTREF(J2).GT.INVDTOACTIVE(CONJ(J2))) INVDTOACTIVE(CONJ(J2))=1.0D0/CONDISTREF(J2)1866:                IF (1.0D0/CONDISTREF(J2).GT.INVDTOACTIVE(CONJ(J2))) INVDTOACTIVE(CONJ(J2))=1.0D0/CONDISTREF(J2)
1877: !              INVDTOACTIVE(CONJ(J2))=INVDTOACTIVE(CONJ(J2))+1.0D0/CONDISTREF(J2)1867: !              INVDTOACTIVE(CONJ(J2))=INVDTOACTIVE(CONJ(J2))+1.0D0/CONDISTREF(J2)
1878:             ENDIF1868:             ENDIF
1879:          ENDIF1869:          ENDIF
1880:          IF (ATOMACTIVE(CONJ(J2))) THEN1870:          IF (ATOMACTIVE(CONJ(J2))) THEN
1881:             IF (.NOT.ATOMACTIVE(CONI(J2))) THEN1871:             IF (.NOT.ATOMACTIVE(CONI(J2))) THEN
1882:                NCONTOACTIVE(CONI(J2))=NCONTOACTIVE(CONI(J2))+11872:                NCONTOACTIVE(CONI(J2))=NCONTOACTIVE(CONI(J2))+1
2176: 543         CONTINUE2166: 543         CONTINUE
2177:          ENDDO2167:          ENDDO
2178:          ATOMACTIVE(NEWATOM)=.TRUE.2168:          ATOMACTIVE(NEWATOM)=.TRUE.
2179:          NACTIVE=NACTIVE+12169:          NACTIVE=NACTIVE+1
2180:          IF (MAXNACTIVE.EQ.0) MAXNACTIVE=NATOMS2170:          IF (MAXNACTIVE.EQ.0) MAXNACTIVE=NATOMS
2181: !2171: !
2182: ! Freeze atoms that became active more than NACTIVE-MAXNACTIVE events ago.2172: ! Freeze atoms that became active more than NACTIVE-MAXNACTIVE events ago.
2183: ! For example, with MAXNACTIVE=5 and 40 active atoms, we would freeze those 2173: ! For example, with MAXNACTIVE=5 and 40 active atoms, we would freeze those 
2184: ! turned on first, second, up to the 35th in the TURNONORDER list.2174: ! turned on first, second, up to the 35th in the TURNONORDER list.
2185: !2175: !
 2176:          IF (DEBUG) WRITE(*,'(A,I6)') 'doaddatom> Number of active atoms is now ',NACTIVE
2186:          IF (NACTIVE.GT.MAXNACTIVE) THEN2177:          IF (NACTIVE.GT.MAXNACTIVE) THEN
2187: !           WRITE(*,'(A)') 'doaddatom> TURNONORDER:'2178: !           WRITE(*,'(A)') 'doaddatom> TURNONORDER:'
2188: !           WRITE(*,'(5I6)') TURNONORDER(1:NACTIVE-1)2179: !           WRITE(*,'(5I6)') TURNONORDER(1:NACTIVE-1)
2189:             NDUMMY=TURNONORDER(NACTIVE-MAXNACTIVE)2180:             NDUMMY=TURNONORDER(NACTIVE-MAXNACTIVE)
2190:             IF (INTFROZEN(NDUMMY)) THEN2181:             IF (INTFROZEN(NDUMMY)) THEN
2191:                IF (DEBUG) WRITE(*,'(A,I6,A,2I6)') ' doaddatom> Not turning off frozen active atom ',NDUMMY,' already frozen'2182:                IF (DEBUG) WRITE(*,'(A,I6,A,2I6)') ' doaddatom> Not turning off frozen active atom ',NDUMMY,' already frozen'
2192:             ELSE2183:             ELSE
2193:                IF (DEBUG) WRITE(*,'(A,I6,A,2I6)') ' doaddatom> Freezing active atom ',NDUMMY2184:                IF (DEBUG) WRITE(*,'(A,I6,A,2I6)') ' doaddatom> Freezing active atom ',NDUMMY
2194:                INTFROZEN(NDUMMY)=.TRUE.2185:                INTFROZEN(NDUMMY)=.TRUE.
2195: !2186: !
2240:             STOP2231:             STOP
2241:          ENDIF2232:          ENDIF
2242: 2233: 
2243:          TURNONORDER(NACTIVE)=NEWATOM2234:          TURNONORDER(NACTIVE)=NEWATOM
2244: !2235: !
2245: ! Initial guess for new active atom position. This is crucial for success in INTCONSTRAINT schemes!2236: ! Initial guess for new active atom position. This is crucial for success in INTCONSTRAINT schemes!
2246: !2237: !
2247:          ESAVED=1.0D1002238:          ESAVED=1.0D100
2248:          ESAVE0=1.0D1002239:          ESAVE0=1.0D100
2249:          ESAVEC=1.0D1002240:          ESAVEC=1.0D100
2250:          FTEST=.TRUE. 
2251:          IF (NCONFORNEWATOM.GE.3) THEN2241:          IF (NCONFORNEWATOM.GE.3) THEN
2252: !2242: !
2253: ! Move the new atom consistently in the local environment of its three nearest actively constrained atoms.2243: ! Move the new atom consistently in the local environment of its three nearest actively constrained atoms.
2254: ! Make a local orthogonal coordinate system and use constant components in this basis.2244: ! Make a local orthogonal coordinate system and use constant components in this basis.
2255: !2245: !
2256:             IF (DEBUG) WRITE(*,'(A,3I6)') ' intlbfgs> initial guess from closest three constrained active atoms, ',CONLIST(1:3)2246:             IF (DEBUG) WRITE(*,'(A,3I6)') ' intlbfgs> initial guess from closest three constrained active atoms, ',CONLIST(1:3)
2257:             VEC1(1:3)=XYZ(3*(CONLIST(2)-1)+1:3*(CONLIST(2)-1)+3)-XYZ(3*(CONLIST(1)-1)+1:3*(CONLIST(1)-1)+3)2247:             VEC1(1:3)=XYZ(3*(CONLIST(2)-1)+1:3*(CONLIST(2)-1)+3)-XYZ(3*(CONLIST(1)-1)+1:3*(CONLIST(1)-1)+3)
2258:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)2248:             DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
2259:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY2249:             IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
2260:             VEC2(1:3)=XYZ(3*(CONLIST(3)-1)+1:3*(CONLIST(3)-1)+3)-XYZ(3*(CONLIST(1)-1)+1:3*(CONLIST(1)-1)+3)2250:             VEC2(1:3)=XYZ(3*(CONLIST(3)-1)+1:3*(CONLIST(3)-1)+3)-XYZ(3*(CONLIST(1)-1)+1:3*(CONLIST(1)-1)+3)
2267:             VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2257:             VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2268:             C1=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC1(1)+ &2258:             C1=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC1(1)+ &
2269:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC1(2)+ &2259:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC1(2)+ &
2270:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC1(3)2260:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC1(3)
2271:             C2=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC2(1)+ &2261:             C2=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC2(1)+ &
2272:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC2(2)+ &2262:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC2(2)+ &
2273:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC2(3)2263:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC2(3)
2274:             C3=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC3(1)+ &2264:             C3=(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))*VEC3(1)+ &
2275:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC3(2)+ &2265:   &            (XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))*VEC3(2)+ &
2276:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC3(3)2266:   &            (XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))*VEC3(3)
2277:  
2278:             DO J1=2,INTIMAGE+12267:             DO J1=2,INTIMAGE+1
2279:                VEC1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3) &2268:                VEC1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3) &
2280:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)2269:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)
2281:                DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)2270:                DUMMY=SQRT(VEC1(1)**2+VEC1(2)**2+VEC1(3)**2)
2282:                IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY2271:                IF (DUMMY.NE.0.0D0) VEC1(1:3)=VEC1(1:3)/DUMMY
2283:                VEC2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3) &2272:                VEC2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3) &
2284:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)2273:   &                     -XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)
2285:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)2274:                DUMMY=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
2286:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)2275:                VEC2(1:3)=VEC2(1:3)-DUMMY*VEC1(1:3)
2287:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)2276:                DUMMY=SQRT(VEC2(1)**2+VEC2(2)**2+VEC2(3)**2)
2288:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY2277:                IF (DUMMY.NE.0.0D0) VEC2(1:3)=VEC2(1:3)/DUMMY
2289:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)2278:                VEC3(1)= VEC1(2)*VEC2(3)-VEC1(3)*VEC2(2)
2290:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)2279:                VEC3(2)=-VEC1(1)*VEC2(3)+VEC1(3)*VEC2(1)
2291:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)2280:                VEC3(3)= VEC1(1)*VEC2(2)-VEC1(2)*VEC2(1)
2292:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &2281:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)= &
2293:   &            XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)+C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)2282:   &            XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)+C1*VEC1(1:3)+C2*VEC2(1:3)+C3*VEC3(1:3)
2294:  
2295: ! 
2296: ! Alternative analytical solution from intersection of three spheres by trilateration 
2297: ! 
2298:                P1(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3) 
2299:                P2(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(2)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(2)-1)+3) 
2300:                P3(1:3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(3)-1)+1:(J1-1)*3*NATOMS+3*(CONLIST(3)-1)+3) 
2301:                R1=CONDIST(1) 
2302:                R2=CONDIST(2) 
2303:                R3=CONDIST(3) 
2304:                CALL TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST) 
2305:                IF (FTEST) THEN 
2306: !                 WRITE(*,'(A,I8)') ' intlbfgs> WARNING *** no trilateration solution for image ',J1 
2307:                ELSE 
2308: !                 WRITE(*,'(A,I8)') ' intlbfgs>                trilateration solution for image ',J1 
2309: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL1=',SOL1(1:3) 
2310: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> SOL2=',SOL2(1:3) 
2311: !                 WRITE(*,'(A,3F20.10)') ' intlbfgs> prev=',XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3) 
2312:                   D1SQ=(SOL1(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2313:   &                   +(SOL1(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2314:   &                   +(SOL1(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2315:                   D2SQ=(SOL2(1)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+1))**2 & 
2316:   &                   +(SOL2(2)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+2))**2 & 
2317:   &                   +(SOL2(3)-XYZ((J1-2)*3*NATOMS+3*(NEWATOM-1)+3))**2 
2318: !                 WRITE(*,'(A,2F20.10)') 'D1SQ,D2SQ=',D1SQ,D2SQ 
2319:                   IF (D1SQ.LT.D2SQ) THEN 
2320:                      XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL1(1:3) 
2321:                   ELSE 
2322:                      XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=SOL2(1:3) 
2323:                   ENDIF 
2324:                ENDIF 
2325:  
2326:             ENDDO2283:             ENDDO
2327:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2284:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2328:             IF (QCIADDREP.GT.0) THEN2285:             IF (QCIADDREP.GT.0) THEN
2329:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2286:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2330:             ELSEIF (CHECKCONINT) THEN2287:             ELSEIF (CHECKCONINT) THEN
2331:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2288:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2332:             ELSE2289:             ELSE
2333:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2290:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2334:             ENDIF2291:             ENDIF
2335:             ESAVE0=ETOTAL2292:             ESAVE0=ETOTAL
2336:             DO J1=2,INTIMAGE+12293:             DO J1=2,INTIMAGE+1
2337:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2294:                XSAVE0(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2338:             ENDDO2295:             ENDDO
2339:          ENDIF2296:          ENDIF
2340:          FTEST=.TRUE.  !!!!!!  DEBUG DJW2297:          IF (NDFORNEWATOM.GE.3) THEN
2341:          IF (FTEST.AND.(NDFORNEWATOM.GE.3)) THEN 
2342: !2298: !
2343: ! Choose three atoms from the BESTPRESERVEDN list at random with bias towards the 2299: ! Choose three atoms from the BESTPRESERVEDN list at random with bias towards the 
2344: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate2300: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate
2345: ! the sum to normalise.2301: ! the sum to normalise.
2346: !2302: !
2347:             DUMMY=0.0D02303:             DUMMY=0.0D0
2348:             DO J1=1,NDFORNEWATOM2304:             DO J1=1,NDFORNEWATOM
2349: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)2305: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)
2350: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTPRESERVEDD(J1))2306: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTPRESERVEDD(J1))
2351:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)2307:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)
2437:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2393:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2438:             ELSE2394:             ELSE
2439:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2395:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2440:             ENDIF2396:             ENDIF
2441:             ESAVED=ETOTAL2397:             ESAVED=ETOTAL
2442:             DO J1=2,INTIMAGE+12398:             DO J1=2,INTIMAGE+1
2443:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2399:                XSAVED(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2444:             ENDDO2400:             ENDDO
2445:          ENDIF2401:          ENDIF
2446: 2402: 
2447:          IF (FTEST.AND.(NCFORNEWATOM.GE.3)) THEN2403:          IF (NCFORNEWATOM.GE.3) THEN
2448: !2404: !
2449: ! Choose three atoms from the BESTCLOSEST list at random with bias towards the2405: ! Choose three atoms from the BESTCLOSEST list at random with bias towards the
2450: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate2406: ! start of the list. Let the relative weight for position i be 1/i**2 and calculate
2451: ! the sum to normalise.2407: ! the sum to normalise.
2452: !2408: !
2453:             DUMMY=0.0D02409:             DUMMY=0.0D0
2454:             DO J1=1,NCFORNEWATOM2410:             DO J1=1,NCFORNEWATOM
2455: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)2411: !              DUMMY=DUMMY+1.0D0/(1.0D0*J1)
2456: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTCLOSESTD(J1))2412: !              DUMMY=DUMMY+1.0D0/(1.0D0*BESTCLOSESTD(J1))
2457:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)2413:                DUMMY=DUMMY+1.0D0/(1.0D0*J1**2)
2532:             DO J1=2,INTIMAGE+12488:             DO J1=2,INTIMAGE+1
2533:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)2489:                XSAVEC(1:3,J1)=XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)
2534:             ENDDO2490:             ENDDO
2535:          ENDIF2491:          ENDIF
2536: !2492: !
2537: ! Standard linear interpolation, with constraint distance scaled by FRAC.2493: ! Standard linear interpolation, with constraint distance scaled by FRAC.
2538: ! Works for FRAC as small as 0.1 with repulsion turned off.2494: ! Works for FRAC as small as 0.1 with repulsion turned off.
2539: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements2495: ! We use an appropriately weighted displacement from atom CONLIST(1) using the displacements
2540: ! in the two end points.2496: ! in the two end points.
2541: !2497: !
2542:          ETOTAL=1.0D62498:          FRAC=1.0D0
2543:          IF (FTEST) THEN2499:          DO J1=2,INTIMAGE+1
2544:             FRAC=1.0D02500:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1)  &
2545:             DO J1=2,INTIMAGE+1 
2546:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+1)  & 
2547:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &2501:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+1)-XYZ(3*(CONLIST(1)-1)+1))/(INTIMAGE+1) &
2548:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)2502:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+1)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+1))/(INTIMAGE+1)
2549:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &2503:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+2)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+2)  &
2550:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &2504:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+2)-XYZ(3*(CONLIST(1)-1)+2))/(INTIMAGE+1) &
2551:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)2505:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+2)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+2))/(INTIMAGE+1)
2552:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &2506:             XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XYZ((J1-1)*3*NATOMS+3*(CONLIST(1)-1)+3)  &
2553:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &2507:  &            +(INTIMAGE-J1+2)*FRAC*(XYZ(3*(NEWATOM-1)+3)-XYZ(3*(CONLIST(1)-1)+3))/(INTIMAGE+1) &
2554:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)2508:  &   +(J1-1)*(XYZ(3*NATOMS*(INTIMAGE+1)+3*(NEWATOM-1)+3)-XYZ(3*NATOMS*(INTIMAGE+1)+3*(CONLIST(1)-1)+3))/(INTIMAGE+1)
2555:             ENDDO2509:          ENDDO
2556:             CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list2510:          CALL CHECKREP(INTIMAGE,XYZ,(3*NATOMS),NNREPSAVE,NREPSAVE+1) ! set up repulsive neighbour list
2557:             IF (QCIADDREP.GT.0) THEN2511:          IF (QCIADDREP.GT.0) THEN
2558:                CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2512:             CALL CONGRAD3(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2559:             ELSEIF (CHECKCONINT) THEN2513:          ELSEIF (CHECKCONINT) THEN
2560:                CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)2514:             CALL CONGRAD2(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2561:             ELSE 
2562:                CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS) 
2563:             ENDIF 
2564:             IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', & 
2565:   &                 ESAVE0,ESAVED,ESAVEC,ETOTAL 
2566:          ELSE2515:          ELSE
2567:             IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> Using trilateration scheme for interpolation'2516:             CALL CONGRAD(NMAXINT,NMININT,ETOTAL,XYZ,GGG,EEE,IMGFREEZE,RMS)
2568:          ENDIF2517:          ENDIF
2569:     2518:          IF (DEBUG) WRITE(*,'(A,4G15.5)') ' intlbfgs> energies for constrained, preserved, closest, and linear schemes=', &
 2519:   &                 ESAVE0,ESAVED,ESAVEC,ETOTAL
2570:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN2520:          IF ((ETOTAL.LT.ESAVEC).AND.(ETOTAL.LT.ESAVED).AND.(ETOTAL.LT.ESAVE0)) THEN
2571:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'2521:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from linear interpolation'
2572:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN2522:          ELSE IF ((ESAVEC.LT.ESAVED).AND.(ESAVEC.LT.ESAVE0)) THEN
2573:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'2523:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest atoms'
2574:             DO J1=2,INTIMAGE+12524:             DO J1=2,INTIMAGE+1
2575:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVEC(1:3,J1)2525:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVEC(1:3,J1)
2576:             ENDDO2526:             ENDDO
2577:             ETOTAL=ESAVEC2527:             ETOTAL=ESAVEC
2578:          ELSE IF (ESAVED.LT.ESAVE0) THEN2528:          ELSE IF (ESAVED.LT.ESAVE0) THEN
2579:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using preserved distances'2529:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using preserved distances'
2580:             DO J1=2,INTIMAGE+12530:             DO J1=2,INTIMAGE+1
2581:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVED(1:3,J1)2531:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVED(1:3,J1)
2582:             ENDDO2532:             ENDDO
2583:             ETOTAL=ESAVED2533:             ETOTAL=ESAVED
2584:          ELSE 2534:          ELSE 
2585:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> interpolation using closest constraints'2535:             IF (DEBUG) WRITE(*,'(A,2G20.10)') ' intlbfgs> lowest energy from interpolation using closest constraints'
2586:             DO J1=2,INTIMAGE+12536:             DO J1=2,INTIMAGE+1
2587:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVE0(1:3,J1)2537:                XYZ((J1-1)*3*NATOMS+3*(NEWATOM-1)+1:(J1-1)*3*NATOMS+3*(NEWATOM-1)+3)=XSAVE0(1:3,J1)
2588:             ENDDO2538:             ENDDO
2589:             ETOTAL=ESAVE02539:             ETOTAL=ESAVE0
2590:          ENDIF2540:          ENDIF
2591:       ENDIF2541:       ENDIF
2592:       NADDED=NADDED+12542:       NADDED=NADDED+1
2593:       IF (NADDED.LT.NTOADD) GOTO 5422543:       IF (NADDED.LT.NTOADD) GOTO 542
2594: !2544: !
2595: ! Check whether we've added all atoms in the amino acid corresponding to the new atom. If not, go back to the top2545: ! Check whether we've added all atoms in the amino acid corresponding to the new atom. If not, go back to the top
3272:       WRITE(*,'(A,2I6,2G15.5,2I6)') ' qcioptperm> images,D12,D21,permutation: ',IM,IM+1,D12,D21,PERM(I1),PERM(I2)3222:       WRITE(*,'(A,2I6,2G15.5,2I6)') ' qcioptperm> images,D12,D21,permutation: ',IM,IM+1,D12,D21,PERM(I1),PERM(I2)
3273:       IF (NSETS(DOGROUP).GT.0) WRITE(*,'(A,6I6)') ' qcioptperm> associated permutation: ', &3223:       IF (NSETS(DOGROUP).GT.0) WRITE(*,'(A,6I6)') ' qcioptperm> associated permutation: ', &
3274:   &                         (PERM(SETS(I1,J3)),PERM(SETS(I2,J3)),J3=1,NSETS(DOGROUP))  3224:   &                         (PERM(SETS(I1,J3)),PERM(SETS(I2,J3)),J3=1,NSETS(DOGROUP))  
3275:    ENDIF3225:    ENDIF
3276: ELSE3226: ELSE
3277:    WRITE(*,'(A,I6,A,I6)') ' qcioptperm> Unknown group size ',NPERMSIZE(DOGROUP),' for group ',DOGROUP3227:    WRITE(*,'(A,I6,A,I6)') ' qcioptperm> Unknown group size ',NPERMSIZE(DOGROUP),' for group ',DOGROUP
3278:    STOP3228:    STOP
3279: ENDIF3229: ENDIF
3280: 3230: 
3281: END SUBROUTINE QCIOPTPERM3231: END SUBROUTINE QCIOPTPERM
3282:  
3283: SUBROUTINE TRILATERATION(P1,P2,P3,R1,R2,R3,SOL1,SOL2,FTEST) 
3284: IMPLICIT NONE 
3285: DOUBLE PRECISION P1(3), P2(3), P3(3), R1, R2, R3, SOL1(3), SOL2(3), I 
3286: DOUBLE PRECISION  TEMP1(3), EX(3), EY(3), EZ(3), DUMMY, TEMP2(3), TEMP3(3), D, J, X, Y, Z, TEMP4 
3287: LOGICAL FTEST 
3288:  
3289: ! # Find the intersection of three spheres                  
3290: ! # P1,P2,P3 are the centers, r1,r2,r3 are the radii        
3291: ! # Implementaton based on Wikipedia Trilateration article.                               
3292:  
3293: FTEST=.FALSE. 
3294: TEMP1(1:3)=P2(1:3)-P1(1:3) 
3295: D=SQRT( TEMP1(1)**2+TEMP1(2)**2+TEMP1(3)**2 ) 
3296: EX(1:3)=TEMP1(1:3)/D 
3297: TEMP2(1:3)=P3(1:3)-P1(1:3) 
3298: I=EX(1)*TEMP2(1)+EX(2)*TEMP2(2)+EX(3)*TEMP2(3) 
3299: TEMP3(1:3)=TEMP2(1:3)-I*EX(1:3) 
3300: DUMMY=SQRT( TEMP3(1)**2+TEMP3(2)**2+TEMP3(3)**2 ) 
3301: EY(1:3)=TEMP3(1:3)/DUMMY 
3302: EZ(1)= EX(2)*EY(3)-EX(3)*EY(2) 
3303: EZ(2)=-EX(1)*EY(3)+EX(3)*EY(1) 
3304: EZ(3)= EX(1)*EY(2)-EX(2)*EY(1) 
3305: J=EY(1)*TEMP2(1)+EY(2)*TEMP2(2)+EY(3)*TEMP2(3) 
3306: X=(R1*R1 - R2*R2 + D*D) / (2.0D0*D) 
3307: Y=(R1*R1 - R3*R3 -2.0D0*I*X + I*I + J*J) / (2.0D0*J) 
3308: TEMP4=R1*R1 - X*X - Y*Y 
3309:  
3310: ! WRITE (*,'(A,9G15.5)') 'trilateration> p1, p2, p3: ',P1(1:3),P2(1:3),P3(1:3) 
3311: ! WRITE (*,'(A,9G15.5)') 'trilateration> ex, ey, ez: ',EX(1:3),EY(1:3),EZ(1:3) 
3312: ! WRITE (*,'(A,9G15.5)') 'trilateration> norms:      ',EX(1)**2+EX(2)**2+EX(3)**2,EY(1)**2+EY(2)**2+EY(3)**2,EZ(1)**2+EZ(2)**2+EZ(3)**2 
3313: ! WRITE (*,'(A,9G15.5)') 'trilateration> r1, r2, r3: ',R1,R2,R3 
3314: ! WRITE (*,'(A,9G15.5)') 'trilateration> X, Y, TEMP4:    ',X,Y,TEMP4 
3315:  
3316: ! PRINT *,'TEMP4=',TEMP4 
3317: !PRINT *,'TEMP4.LT.0.0D0=',TEMP4.LT.0.0D0 
3318:  
3319: IF (TEMP4.LT.0.0D0) THEN 
3320:    FTEST=.TRUE. 
3321:    RETURN 
3322: ELSE 
3323:    FTEST=.FALSE. 
3324:    Z=SQRT(TEMP4) 
3325:    SOL1(1:3)=P1(1:3) + X*EX(1:3) + Y*EY(1:3) + Z*EZ(1:3) 
3326:    SOL2(1:3)=P1(1:3) + X*EX(1:3) + Y*EY(1:3) - Z*EZ(1:3) 
3327: !  PRINT *,'Z=',Z 
3328: ENDIF 
3329:  
3330: END SUBROUTINE TRILATERATION  


r33431/key.f90 2017-11-01 21:30:30.211388375 +0000 r33430/key.f90 2017-11-01 21:30:31.551406207 +0000
 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, QCIRESETINT1, QCIRESETINT2, JMAXCON, NCONOFF 29:      &        MAX_ANGMOM, BNB_NSTEPS, QUIPZ, QCIRESETINT1, QCIRESETINT2, JMAXCON
 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,&
180:       DOUBLE PRECISION, ALLOCATABLE :: LJADDREP(:,:), LJADDATT(:,:)180:       DOUBLE PRECISION, ALLOCATABLE :: LJADDREP(:,:), LJADDATT(:,:)
181:       LOGICAL, ALLOCATABLE :: CONACTIVE(:)181:       LOGICAL, ALLOCATABLE :: CONACTIVE(:)
182:       LOGICAL, ALLOCATABLE :: ATOMACTIVE(:)182:       LOGICAL, ALLOCATABLE :: ATOMACTIVE(:)
183:       INTEGER, ALLOCATABLE :: NITSTART(:)183:       INTEGER, ALLOCATABLE :: NITSTART(:)
184:       DOUBLE PRECISION, ALLOCATABLE :: RPMASSES(:), XMINA(:), XMINB(:)184:       DOUBLE PRECISION, ALLOCATABLE :: RPMASSES(:), XMINA(:), XMINB(:)
185:       DOUBLE PRECISION, ALLOCATABLE :: RBOPS(:,:)185:       DOUBLE PRECISION, ALLOCATABLE :: RBOPS(:,:)
186:       DOUBLE PRECISION, ALLOCATABLE :: SAVES(:), SAVEF(:)186:       DOUBLE PRECISION, ALLOCATABLE :: SAVES(:), SAVEF(:)
187: !     LOGICAL, ALLOCATABLE :: CONTEST(:,:)187: !     LOGICAL, ALLOCATABLE :: CONTEST(:,:)
188:       INTEGER, ALLOCATABLE :: ORDERI(:), ORDERJ(:), REPPOW(:)188:       INTEGER, ALLOCATABLE :: ORDERI(:), ORDERJ(:), REPPOW(:)
189:       INTEGER, ALLOCATABLE :: CONI(:), CONJ(:), CONION(:), CONJON(:), CONOFFLIST(:)189:       INTEGER, ALLOCATABLE :: CONI(:), CONJ(:), CONION(:), CONJON(:), CONOFFLIST(:)
190:       LOGICAL, ALLOCATABLE :: CONOFFTRIED(:) 
191:       INTEGER, ALLOCATABLE :: CONIFIX(:), CONJFIX(:), REPIFIX(:), REPJFIX(:)190:       INTEGER, ALLOCATABLE :: CONIFIX(:), CONJFIX(:), REPIFIX(:), REPJFIX(:)
192:       INTEGER, ALLOCATABLE :: REPI(:), REPJ(:)191:       INTEGER, ALLOCATABLE :: REPI(:), REPJ(:)
193:       INTEGER, ALLOCATABLE :: CPCONI(:), CPCONJ(:)192:       INTEGER, ALLOCATABLE :: CPCONI(:), CPCONJ(:)
194:       INTEGER, ALLOCATABLE :: CPREPI(:), CPREPJ(:)193:       INTEGER, ALLOCATABLE :: CPREPI(:), CPREPJ(:)
195:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)194:       DOUBLE PRECISION, ALLOCATABLE :: REPCUT(:), NREPCUT(:), CPREPCUT(:), REPCUTFIX(:)
196:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)195:       INTEGER, ALLOCATABLE :: NREPI(:), NREPJ(:)
197:       INTEGER, ALLOCATABLE :: BESTPERM(:)196:       INTEGER, ALLOCATABLE :: BESTPERM(:)
198:       INTEGER, ALLOCATABLE :: RBGROUP(:), RBNINGROUP(:)197:       INTEGER, ALLOCATABLE :: RBGROUP(:), RBNINGROUP(:)
199:       INTEGER, ALLOCATABLE :: TO_ALIGN(:)198:       INTEGER, ALLOCATABLE :: TO_ALIGN(:)
200:       LOGICAL :: SETCHIRAL=.FALSE.199:       LOGICAL :: SETCHIRAL=.FALSE.


r33431/keywords.f 2017-11-01 21:30:30.439391407 +0000 r33430/keywords.f 2017-11-01 21:30:31.839410037 +0000
3619:                READ(LUNIT,*) CONCUTFIX(1:NCONSTRAINTFIX)3619:                READ(LUNIT,*) CONCUTFIX(1:NCONSTRAINTFIX)
3620:                READ(LUNIT,*) NREPULSIVEFIX3620:                READ(LUNIT,*) NREPULSIVEFIX
3621:                ALLOCATE(REPIFIX(NREPULSIVEFIX),REPJFIX(NREPULSIVEFIX),REPCUTFIX(NREPULSIVEFIX))3621:                ALLOCATE(REPIFIX(NREPULSIVEFIX),REPJFIX(NREPULSIVEFIX),REPCUTFIX(NREPULSIVEFIX))
3622:                READ(LUNIT,*) REPIFIX(1:NREPULSIVEFIX)3622:                READ(LUNIT,*) REPIFIX(1:NREPULSIVEFIX)
3623:                READ(LUNIT,*) REPJFIX(1:NREPULSIVEFIX)3623:                READ(LUNIT,*) REPJFIX(1:NREPULSIVEFIX)
3624:                READ(LUNIT,*) REPCUTFIX(1:NREPULSIVEFIX)3624:                READ(LUNIT,*) REPCUTFIX(1:NREPULSIVEFIX)
3625:                CLOSE(LUNIT)3625:                CLOSE(LUNIT)
3626:                PRINT '(A)',' keyword> Constraint potential parameters read from file congeom.dat'3626:                PRINT '(A)',' keyword> Constraint potential parameters read from file congeom.dat'
3627:                INTCONMAX=NCONSTRAINTFIX3627:                INTCONMAX=NCONSTRAINTFIX
3628:                NREPMAX=NREPULSIVEFIX3628:                NREPMAX=NREPULSIVEFIX
3629:                ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),CONOFFTRIED(INTCONMAX))3629:                ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX))
3630:                CONOFFTRIED(1:INTCONMAX)=.FALSE. 
3631:                ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))3630:                ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
3632:                ALLOCATE(CONACTIVE(NCONSTRAINTFIX))3631:                ALLOCATE(CONACTIVE(NCONSTRAINTFIX))
3633:             ELSE3632:             ELSE
3634:                INQUIRE(FILE='congeom',EXIST=CONFILE)3633:                INQUIRE(FILE='congeom',EXIST=CONFILE)
3635:                NCONGEOM=03634:                NCONGEOM=0
3636:                IF (.NOT.CONFILE) THEN3635:                IF (.NOT.CONFILE) THEN
3637:                   PRINT '(A)',' keyword> WARNING *** no congeom file found. Will use end point minima only.'3636:                   PRINT '(A)',' keyword> WARNING *** no congeom file found. Will use end point minima only.'
3638:                ELSE3637:                ELSE
3639:                   LUNIT=GETUNIT()3638:                   LUNIT=GETUNIT()
3640:                   OPEN(LUNIT,FILE='congeom',STATUS='OLD')3639:                   OPEN(LUNIT,FILE='congeom',STATUS='OLD')


r33431/make_conpot.f90 2017-11-01 21:30:30.659394335 +0000 r33430/make_conpot.f90 2017-11-01 21:30:32.059412966 +0000
 11: !   GNU General Public License for more details. 11: !   GNU General Public License for more details.
 12: ! 12: !
 13: !   You should have received a copy of the GNU General Public License 13: !   You should have received a copy of the GNU General Public License
 14: !   along with this program; if not, write to the Free Software 14: !   along with this program; if not, write to the Free Software
 15: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 15: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 16: ! 16: !
 17: SUBROUTINE MAKE_CONPOT(NCPFIT,MINCOORDS) 17: SUBROUTINE MAKE_CONPOT(NCPFIT,MINCOORDS)
 18: USE KEY, ONLY : INTCONSEP, NREPMAX, NREPULSIVE, CONDISTREF, REPCON, INTCONSTRAINTREP, & 18: USE KEY, ONLY : INTCONSEP, NREPMAX, NREPULSIVE, CONDISTREF, REPCON, INTCONSTRAINTREP, &
 19:   & REPCUT, NCONSTRAINT, CONI, CONJ, CONDISTREFLOCAL, INTCONMAX, CONACTIVE, & 19:   & REPCUT, NCONSTRAINT, CONI, CONJ, CONDISTREFLOCAL, INTCONMAX, CONACTIVE, &
 20:   & INTCONSTRAINREPCUT, INTREPSEP, REPI, REPJ, INTCONSTRAINTTOL, REPCUT, NREPI, NREPJ, NREPCUT, & 20:   & INTCONSTRAINREPCUT, INTREPSEP, REPI, REPJ, INTCONSTRAINTTOL, REPCUT, NREPI, NREPJ, NREPCUT, &
 21:   & NCONGEOM, CONGEOM, NNREPULSIVE, BULKT, RIGIDBODY, TWOD, CONOFFLIST, CONOFFTRIED, & 21:   & NCONGEOM, CONGEOM, NNREPULSIVE, BULKT, RIGIDBODY, TWOD, CONOFFLIST, &
 22:   & INTFROZEN, FREEZE, INTFREEZET, INTFREEZETOL, INTFREEZEMIN, CONIFIX, CONJFIX, CONDISTREFFIX, REPIFIX, REPJFIX, & 22:   & INTFROZEN, FREEZE, INTFREEZET, INTFREEZETOL, INTFREEZEMIN, CONIFIX, CONJFIX, CONDISTREFFIX, REPIFIX, REPJFIX, &
 23:   & REPCUTFIX, NCONGEOM, NREPULSIVEFIX, CONDATT, NCONSTRAINTFIX, CONCUTLOCAL, CONCUTFIX, CONCUT 23:   & REPCUTFIX, NCONGEOM, NREPULSIVEFIX, CONDATT, NCONSTRAINTFIX, CONCUTLOCAL, CONCUTFIX, CONCUT
 24: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3 24: USE COMMONS, ONLY: NATOMS, DEBUG, PARAM1, PARAM2, PARAM3
 25:  25: 
 26: IMPLICIT NONE  26: IMPLICIT NONE 
 27: DOUBLE PRECISION DF, D, RMAT(3,3), DISTANCE, D2 27: DOUBLE PRECISION DF, D, RMAT(3,3), DISTANCE, D2
 28: INTEGER :: J2,ISTAT,J1,J3,J4,NCPFIT,J5,NQCIFREEZE,NDUMMY,LUNIT,GETUNIT 28: INTEGER :: J2,ISTAT,J1,J3,J4,NCPFIT,J5,NQCIFREEZE,NDUMMY,LUNIT,GETUNIT
 29: INTEGER NCONFORNEWATOM 29: INTEGER NCONFORNEWATOM
 30: DOUBLE PRECISION :: NDIST, MINCOORDS(NCPFIT,3*NATOMS), DMIN, LINTCONSTRAINTTOL, & 30: DOUBLE PRECISION :: NDIST, MINCOORDS(NCPFIT,3*NATOMS), DMIN, LINTCONSTRAINTTOL, &
 31:   &                 LXYZ(6*NATOMS) 31:   &                 LXYZ(6*NATOMS)
102: 102: 
103: IF (NATOMS-NQCIFREEZE.LT.INTFREEZEMIN) THEN103: IF (NATOMS-NQCIFREEZE.LT.INTFREEZEMIN) THEN
104:    DO J1=NATOMS,NATOMS-INTFREEZEMIN+1,-1104:    DO J1=NATOMS,NATOMS-INTFREEZEMIN+1,-1
105:       INTFROZEN(DLIST(J1))=.FALSE.105:       INTFROZEN(DLIST(J1))=.FALSE.
106:    ENDDO106:    ENDDO
107:    NQCIFREEZE=MAX(0,NATOMS-INTFREEZEMIN)107:    NQCIFREEZE=MAX(0,NATOMS-INTFREEZEMIN)
108:    IF (DEBUG) WRITE(*, '(A,I6,A)') ' make_conpot> Freezing ',NQCIFREEZE,' atoms'108:    IF (DEBUG) WRITE(*, '(A,I6,A)') ' make_conpot> Freezing ',NQCIFREEZE,' atoms'
109: ENDIF109: ENDIF
110: 110: 
111: IF (.NOT.ALLOCATED(CONI)) THEN 111: IF (.NOT.ALLOCATED(CONI)) THEN 
112:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),CONOFFTRIED(INTCONMAX))112:    ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX))
113:    CONOFFTRIED(1:INTCONMAX)=.FALSE. 
114:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))113:    ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
115: ENDIF114: ENDIF
116: 115: 
117: IF (NQCIFREEZE.EQ.NATOMS) THEN116: IF (NQCIFREEZE.EQ.NATOMS) THEN
118:    NREPULSIVE=0117:    NREPULSIVE=0
119:    NNREPULSIVE=0118:    NNREPULSIVE=0
120:    NCONSTRAINT=0119:    NCONSTRAINT=0
121:    IF (DEBUG) WRITE(*,'(A,2I10,A,G20.10)') ' make_conpot> Total number of constraints and repulsions=', &120:    IF (DEBUG) WRITE(*,'(A,2I10,A,G20.10)') ' make_conpot> Total number of constraints and repulsions=', &
122:   &   NCONSTRAINT,NREPULSIVE121:   &   NCONSTRAINT,NREPULSIVE
123:    122:    
469: NREPCUT(1:NREPMAX)=REPTEMP(1:NREPMAX)468: NREPCUT(1:NREPMAX)=REPTEMP(1:NREPMAX)
470: 469: 
471: DEALLOCATE(IREPTEMP,REPTEMP)470: DEALLOCATE(IREPTEMP,REPTEMP)
472: NREPMAX=2*NREPMAX471: NREPMAX=2*NREPMAX
473: 472: 
474: END SUBROUTINE REPDOUBLE473: END SUBROUTINE REPDOUBLE
475: 474: 
476: 475: 
477: SUBROUTINE CONDOUBLE476: SUBROUTINE CONDOUBLE
478: USE KEY, ONLY : CONI, CONJ, CONDISTREF, INTCONMAX, CONIFIX, CONJFIX, CONDISTREFFIX, NCONGEOM, &477: USE KEY, ONLY : CONI, CONJ, CONDISTREF, INTCONMAX, CONIFIX, CONJFIX, CONDISTREFFIX, NCONGEOM, &
479:   &             CONCUT, CONCUTFIX, CONOFFLIST, CONOFFTRIED478:   &             CONCUT, CONCUTFIX, CONOFFLIST
480: IMPLICIT NONE479: IMPLICIT NONE
481: DOUBLE PRECISION, ALLOCATABLE :: CPTEMP(:)480: DOUBLE PRECISION, ALLOCATABLE :: CPTEMP(:)
482: INTEGER, ALLOCATABLE :: ICPTEMP(:)481: INTEGER, ALLOCATABLE :: ICPTEMP(:)
483: LOGICAL, ALLOCATABLE :: LCPTEMP(:) 
484: 482: 
485: ALLOCATE(ICPTEMP(INTCONMAX))483: ALLOCATE(ICPTEMP(INTCONMAX))
486: ALLOCATE(CPTEMP(INTCONMAX))484: ALLOCATE(CPTEMP(1:INTCONMAX))
487: ALLOCATE(LCPTEMP(INTCONMAX)) 
488:                 485:                 
489: ICPTEMP(1:INTCONMAX)=CONI(1:INTCONMAX)486: ICPTEMP(1:INTCONMAX)=CONI(1:INTCONMAX)
490: DEALLOCATE(CONI)487: DEALLOCATE(CONI)
491: ALLOCATE(CONI(2*INTCONMAX))488: ALLOCATE(CONI(2*INTCONMAX))
492: CONI(1:INTCONMAX)=ICPTEMP(1:INTCONMAX)489: CONI(1:INTCONMAX)=ICPTEMP(1:INTCONMAX)
493: 490: 
494: IF (NCONGEOM.GE.2) THEN491: IF (NCONGEOM.GE.2) THEN
495:    ICPTEMP(1:INTCONMAX)=CONIFIX(1:INTCONMAX)492:    ICPTEMP(1:INTCONMAX)=CONIFIX(1:INTCONMAX)
496:    DEALLOCATE(CONIFIX)493:    DEALLOCATE(CONIFIX)
497:    ALLOCATE(CONIFIX(2*INTCONMAX))494:    ALLOCATE(CONIFIX(2*INTCONMAX))
516: ICPTEMP(1:INTCONMAX)=CONJ(1:INTCONMAX)513: ICPTEMP(1:INTCONMAX)=CONJ(1:INTCONMAX)
517: DEALLOCATE(CONJ)514: DEALLOCATE(CONJ)
518: ALLOCATE(CONJ(2*INTCONMAX))515: ALLOCATE(CONJ(2*INTCONMAX))
519: CONJ(1:INTCONMAX)=ICPTEMP(1:INTCONMAX)516: CONJ(1:INTCONMAX)=ICPTEMP(1:INTCONMAX)
520:                517:                
521: ICPTEMP(1:INTCONMAX)=CONOFFLIST(1:INTCONMAX)518: ICPTEMP(1:INTCONMAX)=CONOFFLIST(1:INTCONMAX)
522: DEALLOCATE(CONOFFLIST)519: DEALLOCATE(CONOFFLIST)
523: ALLOCATE(CONOFFLIST(2*INTCONMAX))520: ALLOCATE(CONOFFLIST(2*INTCONMAX))
524: CONOFFLIST(1:INTCONMAX)=ICPTEMP(1:INTCONMAX)521: CONOFFLIST(1:INTCONMAX)=ICPTEMP(1:INTCONMAX)
525:                522:                
526: LCPTEMP(1:INTCONMAX)=CONOFFTRIED(1:INTCONMAX) 
527: DEALLOCATE(CONOFFTRIED) 
528: ALLOCATE(CONOFFTRIED(2*INTCONMAX)) 
529: CONOFFTRIED(1:INTCONMAX)=LCPTEMP(1:INTCONMAX) 
530: CONOFFTRIED(INTCONMAX+1:2*INTCONMAX)=.FALSE. 
531:                 
532: CPTEMP(1:INTCONMAX)=CONDISTREF(1:INTCONMAX)523: CPTEMP(1:INTCONMAX)=CONDISTREF(1:INTCONMAX)
533: DEALLOCATE(CONDISTREF)524: DEALLOCATE(CONDISTREF)
534: ALLOCATE(CONDISTREF(2*INTCONMAX))525: ALLOCATE(CONDISTREF(2*INTCONMAX))
535: CONDISTREF(1:INTCONMAX)=CPTEMP(1:INTCONMAX)526: CONDISTREF(1:INTCONMAX)=CPTEMP(1:INTCONMAX)
536:                527:                
537: CPTEMP(1:INTCONMAX)=CONCUT(1:INTCONMAX)528: CPTEMP(1:INTCONMAX)=CONCUT(1:INTCONMAX)
538: DEALLOCATE(CONCUT)529: DEALLOCATE(CONCUT)
539: ALLOCATE(CONCUT(2*INTCONMAX))530: ALLOCATE(CONCUT(2*INTCONMAX))
540: CONCUT(1:INTCONMAX)=CPTEMP(1:INTCONMAX)531: CONCUT(1:INTCONMAX)=CPTEMP(1:INTCONMAX)
541: 532: 
542: INTCONMAX=2*INTCONMAX533: INTCONMAX=2*INTCONMAX
543: DEALLOCATE(CPTEMP)534: DEALLOCATE(CPTEMP)
544: DEALLOCATE(ICPTEMP)535: DEALLOCATE(ICPTEMP)
545: DEALLOCATE(LCPTEMP) 
546: 536: 
547: END SUBROUTINE CONDOUBLE537: END SUBROUTINE CONDOUBLE


r33431/OPTIM.F 2017-11-01 21:30:29.543379441 +0000 r33430/OPTIM.F 2017-11-01 21:30:30.883397316 +0000
167: 167: 
168:       IF (AMBERT .OR. NABT) CALL SET_CHECK_CHIRAL(Q,NATOMS,GOODSTRUCTURE1,CHIARRAY1)168:       IF (AMBERT .OR. NABT) CALL SET_CHECK_CHIRAL(Q,NATOMS,GOODSTRUCTURE1,CHIARRAY1)
169: 169: 
170:       IF (CHECKDT) CALL CHECKDRVTS(Q)170:       IF (CHECKDT) CALL CHECKDRVTS(Q)
171: 171: 
172:       IF ((INTCONSTRAINTT.OR.INTLJT).AND.(NCONGEOM.GE.2)) THEN172:       IF ((INTCONSTRAINTT.OR.INTLJT).AND.(NCONGEOM.GE.2)) THEN
173: !173: !
174: ! Set up all the constraints and repulsions for zero frozen atoms.174: ! Set up all the constraints and repulsions for zero frozen atoms.
175: !175: !
176:          IF (.NOT.ALLOCATED(CONI)) THEN176:          IF (.NOT.ALLOCATED(CONI)) THEN
177:             ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX),CONOFFTRIED(INTCONMAX))177:             ALLOCATE(CONI(INTCONMAX),CONJ(INTCONMAX),CONDISTREF(INTCONMAX),CONCUT(INTCONMAX),CONOFFLIST(INTCONMAX))
178:             CONOFFTRIED(1:INTCONMAX)=.FALSE. 
179:             ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))178:             ALLOCATE(REPI(NREPMAX),REPJ(NREPMAX),NREPI(NREPMAX),NREPJ(NREPMAX),REPCUT(NREPMAX),NREPCUT(NREPMAX))
180:          ENDIF179:          ENDIF
181:          INTFREEZETOLSAVE=INTFREEZETOL180:          INTFREEZETOLSAVE=INTFREEZETOL
182:          INTFREEZETOL=-1.0D0181:          INTFREEZETOL=-1.0D0
183:          PRINT *,'OPTIM> NCONGEOM=',NCONGEOM182:          PRINT *,'OPTIM> NCONGEOM=',NCONGEOM
184:          CALL MAKE_CONPOT(NCONGEOM,CONGEOM) 183:          CALL MAKE_CONPOT(NCONGEOM,CONGEOM) 
185:          INTFREEZETOL=INTFREEZETOLSAVE184:          INTFREEZETOL=INTFREEZETOLSAVE
186:       ENDIF185:       ENDIF
187: 186: 
188: 187: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0