hdiff output

r32722/commons.f90 2017-06-05 15:30:19.963789693 +0100 r32721/commons.f90 2017-06-05 15:30:21.955815397 +0100
650:       INTEGER, ALLOCATABLE ::  MLPOUTCOME(:)650:       INTEGER, ALLOCATABLE ::  MLPOUTCOME(:)
651:       DOUBLE PRECISION, ALLOCATABLE ::  MLQDAT(:,:)651:       DOUBLE PRECISION, ALLOCATABLE ::  MLQDAT(:,:)
652:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)652:       INTEGER, ALLOCATABLE ::  MLQOUTCOME(:)
653:       INTEGER, ALLOCATABLE ::  LJADDNN(:,:)653:       INTEGER, ALLOCATABLE ::  LJADDNN(:,:)
654: 654: 
655:       INTEGER, DIMENSION(:,:), ALLOCATABLE :: BONDS !for QCIAMBER655:       INTEGER, DIMENSION(:,:), ALLOCATABLE :: BONDS !for QCIAMBER
656: 656: 
657: !OPEP interface657: !OPEP interface
658:       LOGICAL :: OPEPT, OPEP_RNAT658:       LOGICAL :: OPEPT, OPEP_RNAT
659: 659: 
660: !Orbital variables 
661:       LOGICAL :: ORBITALS 
662:       INTEGER :: NROTS, NORBS, ORBVAREXPONENT 
663:       DOUBLE PRECISION, ALLOCATABLE :: R2INTS(:,:), DIPINTS(:,:,:) 
664: END MODULE COMMONS660: END MODULE COMMONS


r32722/finalio.f90 2017-06-05 15:30:20.183792533 +0100 r32721/finalio.f90 2017-06-05 15:30:22.179818287 +0100
767:         ENDDO767:         ENDDO
768: 768: 
769:      ELSE IF (GBT.OR.GBDT.OR.GBDPT.OR.MSGBT) THEN769:      ELSE IF (GBT.OR.GBDT.OR.GBDPT.OR.MSGBT) THEN
770:         DO J2 = 1, NATOMS/2770:         DO J2 = 1, NATOMS/2
771:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)771:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
772:         ENDDO772:         ENDDO
773:         DO J2 = 1, NATOMS/2773:         DO J2 = 1, NATOMS/2
774:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)774:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*NATOMS/2+3*(J2-1)+J3),J3=1,3)
775:         ENDDO775:         ENDDO
776: 776: 
777:      ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT.OR.MLPVB3T.OR.ORBITALS) THEN777:      ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT.OR.MLPVB3T) THEN
778:         DO J2 = 1, NATOMS778:         DO J2 = 1, NATOMS
779:            WRITE(MYUNIT2,'(3G20.10)') QMINP(J1,J2)779:            WRITE(MYUNIT2,'(3G20.10)') QMINP(J1,J2)
780:         ENDDO780:         ENDDO
781: 781: 
782:      ELSE IF (GEMT) THEN782:      ELSE IF (GEMT) THEN
783:         DO J2 = 1, NATOMS783:         DO J2 = 1, NATOMS
784:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)784:            WRITE(MYUNIT2,'(3f20.10)') (QMINP(J1,3*(J2-1)+J3),J3=1,3)
785:         ENDDO785:         ENDDO
786: 786: 
787:      ELSE IF (BLNT.AND.(.NOT.P46).AND.(.NOT.G46)) THEN787:      ELSE IF (BLNT.AND.(.NOT.P46).AND.(.NOT.G46)) THEN


r32722/io1.f 2017-06-05 15:30:20.403795370 +0100 r32721/io1.f 2017-06-05 15:30:22.499822416 +0100
115:                    READ(COORDS_UNIT,*) COORDS(J1,JP)115:                    READ(COORDS_UNIT,*) COORDS(J1,JP)
116:                ENDDO116:                ENDDO
117:             ENDDO117:             ENDDO
118:          ELSEIF (MLQT) THEN ! for this ML it is one variable per line118:          ELSEIF (MLQT) THEN ! for this ML it is one variable per line
119:             REWIND(COORDS_UNIT)119:             REWIND(COORDS_UNIT)
120:             DO JP=1,NPAR120:             DO JP=1,NPAR
121:                DO J1=1,NATOMS121:                DO J1=1,NATOMS
122:                    READ(COORDS_UNIT,*) COORDS(J1,JP)122:                    READ(COORDS_UNIT,*) COORDS(J1,JP)
123:                ENDDO123:                ENDDO
124:             ENDDO124:             ENDDO
125:          ELSEIF (ORBITALS) THEN ! for this orbital landscape it is one variable per line 
126:             REWIND(COORDS_UNIT) 
127:             DO JP=1,NPAR 
128:                DO J1=1,NATOMS 
129:                    READ(COORDS_UNIT,*) COORDS(J1,JP) 
130:                ENDDO 
131:             ENDDO 
132:          ELSE125:          ELSE
133:             REWIND(COORDS_UNIT)126:             REWIND(COORDS_UNIT)
134:             DO JP=1,NPAR127:             DO JP=1,NPAR
135:                DO J1=1,NATOMS128:                DO J1=1,NATOMS
136:                   J2=3*(J1-1)129:                   J2=3*(J1-1)
137:                    READ(COORDS_UNIT,*) COORDS(J2+1,JP), COORDS(J2+2,JP), COORDS(J2+3,JP)130:                    READ(COORDS_UNIT,*) COORDS(J2+1,JP), COORDS(J2+2,JP), COORDS(J2+3,JP)
138:                ENDDO131:                ENDDO
139:             ENDDO132:             ENDDO
140:          ENDIF133:          ENDIF
141:          CLOSE(COORDS_UNIT)134:          CLOSE(COORDS_UNIT)
202: 195: 
203:       IF (.NOT.SEEDT.AND..NOT.AMHT.AND..NOT.SUPPRESST) THEN196:       IF (.NOT.SEEDT.AND..NOT.AMHT.AND..NOT.SUPPRESST) THEN
204:          WRITE(MYUNIT,20) 197:          WRITE(MYUNIT,20) 
205: 20       FORMAT('Initial coordinates:')198: 20       FORMAT('Initial coordinates:')
206:          IF (MPIT) THEN199:          IF (MPIT) THEN
207:             WRITE(MYUNIT,30) (COORDS(J1,MYNODE+1),J1=1,3*NATOMS)200:             WRITE(MYUNIT,30) (COORDS(J1,MYNODE+1),J1=1,3*NATOMS)
208:          ELSEIF (MLP3T.OR.MLPVB3T) THEN 201:          ELSEIF (MLP3T.OR.MLPVB3T) THEN 
209:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)202:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)
210:          ELSEIF (MLQT) THEN 203:          ELSEIF (MLQT) THEN 
211:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)204:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS)
212:          ELSEIF (ORBITALS) THEN  
213:             WRITE(MYUNIT,'(G20.10)') (COORDS(J1,MYNODE+1),J1=1,NATOMS) 
214:          ELSE 205:          ELSE 
215:            DO JP=1,NPAR206:            DO JP=1,NPAR
216:                WRITE(MYUNIT,30) (COORDS(J1,JP),J1=1,3*NATOMS)207:                WRITE(MYUNIT,30) (COORDS(J1,JP),J1=1,3*NATOMS)
217: 30             FORMAT(3F20.10)208: 30             FORMAT(3F20.10)
218:             ENDDO209:             ENDDO
219:          ENDIF210:          ENDIF
220:       ENDIF211:       ENDIF
221: 212: 
222:       IF (MSORIGT) THEN213:       IF (MSORIGT) THEN
223:          WRITE(MYUNIT,'(I4,A)') NATOMS,' M and S silicon atoms'214:          WRITE(MYUNIT,'(I4,A)') NATOMS,' M and S silicon atoms'
728:             WRITE(MYUNIT,'(2(I4,A))') NRIGIDBODY-NHEXAMERS,' pentamers and ',NHEXAMERS,' hexamers'719:             WRITE(MYUNIT,'(2(I4,A))') NRIGIDBODY-NHEXAMERS,' pentamers and ',NHEXAMERS,' hexamers'
729:             WRITE(MYUNIT,'(A,4G20.10)') 'rho, eps, sigma and pentamer radius=',CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT720:             WRITE(MYUNIT,'(A,4G20.10)') 'rho, eps, sigma and pentamer radius=',CAPSIDRHO,CAPSIDEPS,SIGMAPENT,RADPENT
730:             WRITE(MYUNIT,'(A,4G20.10)') 'hexamer sigma, radius and hex/pent sigma=',SIGMAHEX,RADHEX,SIGMAPH721:             WRITE(MYUNIT,'(A,4G20.10)') 'hexamer sigma, radius and hex/pent sigma=',SIGMAHEX,RADHEX,SIGMAPH
731:          ENDIF722:          ENDIF
732:       ELSEIF (MLPVB3T) THEN723:       ELSEIF (MLPVB3T) THEN
733:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLPVB3'724:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLPVB3'
734:       ELSEIF (MLP3T) THEN725:       ELSEIF (MLP3T) THEN
735:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLP3'726:          WRITE(MYUNIT,'(I4,A)') NATOMS,' link weights for MLP3'
736:       ELSEIF (MLQT) THEN727:       ELSEIF (MLQT) THEN
737:          WRITE(MYUNIT,'(I4,A)') NATOMS,' variables for ML quadratic'728:          WRITE(MYUNIT,'(I4,A)') NATOMS,' variables for ML quadratic'
738:       ELSEIF (ORBITALS) THEN 
739:          WRITE(MYUNIT,'(I4,A)') NATOMS,' rotations for orbital landscape' 
740:       ELSEIF  (LJADDT) THEN729:       ELSEIF  (LJADDT) THEN
741:          IF (SORTT) THEN730:          IF (SORTT) THEN
742:             WRITE(MYUNIT,'(A)') 'Turning off SORT option for LJADD'731:             WRITE(MYUNIT,'(A)') 'Turning off SORT option for LJADD'
743:             SORTT=.FALSE.732:             SORTT=.FALSE.
744:          ENDIF733:          ENDIF
745:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ addressable atoms'734:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ addressable atoms'
746:       ELSE735:       ELSE
747:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ atoms'736:          WRITE(MYUNIT,'(I4,A)') NATOMS,' LJ atoms'
748:       ENDIF737:       ENDIF
749:       IF (PYGPERIODICT.OR.PYBINARYT) CALL INITIALISEPYGPERIODIC738:       IF (PYGPERIODICT.OR.PYBINARYT) CALL INITIALISEPYGPERIODIC
764:          WRITE(MYUNIT, '(A,I6,A)') 'Searching for approximate symmetry elements every ',NSYMINTERVAL,' steps'753:          WRITE(MYUNIT, '(A,I6,A)') 'Searching for approximate symmetry elements every ',NSYMINTERVAL,' steps'
765:          WRITE(MYUNIT, '(A,5F15.5)') 'Distance tolerances: ',SYMTOL1,SYMTOL2,SYMTOL3,SYMTOL4,SYMTOL5754:          WRITE(MYUNIT, '(A,5F15.5)') 'Distance tolerances: ',SYMTOL1,SYMTOL2,SYMTOL3,SYMTOL4,SYMTOL5
766:          WRITE(MYUNIT, '(A,F15.5)') 'Threshold for distinguishing transformation matrices: ',MATDIFF755:          WRITE(MYUNIT, '(A,F15.5)') 'Threshold for distinguishing transformation matrices: ',MATDIFF
767:          WRITE(MYUNIT,'(A,F15.5)') 'Exponential factor in core-weighted centre of mass calculation: ',DISTFAC756:          WRITE(MYUNIT,'(A,F15.5)') 'Exponential factor in core-weighted centre of mass calculation: ',DISTFAC
768:          WRITE(MYUNIT, '(A,I5)') 'Maximum number of quenches for floater/hole permutations=',NSYMQMAX757:          WRITE(MYUNIT, '(A,I5)') 'Maximum number of quenches for floater/hole permutations=',NSYMQMAX
769:          IF (MOVESHELLT) WRITE(MYUNIT,'(A,I8,A,F12.5)') 'Shell moves allowed in blocks of ',SHELLMOVEMAX,' with probability ',758:          IF (MOVESHELLT) WRITE(MYUNIT,'(A,I8,A,F12.5)') 'Shell moves allowed in blocks of ',SHELLMOVEMAX,' with probability ',
770:      &                        SHELLPROB759:      &                        SHELLPROB
771:       ENDIF760:       ENDIF
772:       IF (DEBUG.OR.CHECKMARKOVT) WRITE(MYUNIT,'(A,I6,A)') 'io1> checking the energy of the saved coordinates in the chain'761:       IF (DEBUG.OR.CHECKMARKOVT) WRITE(MYUNIT,'(A,I6,A)') 'io1> checking the energy of the saved coordinates in the chain'
773:       IF (FREEZE) THEN762:       IF (FREEZE) THEN
774:          IF (MLQT.OR.MLPVB3T.OR.MLP3T.OR.ORBITALS) THEN763:          IF (MLQT.OR.MLPVB3T.OR.MLP3T) THEN
775:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' variables will be frozen:'764:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' variables will be frozen:'
776:          ELSE765:          ELSE
777:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' atoms will be frozen:'766:             WRITE(MYUNIT,'(A,I6,A)') 'io1> ', NFREEZE,' atoms will be frozen:'
778:          ENDIF767:          ENDIF
779:          DO J1=1,NATOMS768:          DO J1=1,NATOMS
780:             IF (FROZEN(J1)) WRITE(MYUNIT,'(I6)') J1769:             IF (FROZEN(J1)) WRITE(MYUNIT,'(I6)') J1
781:          ENDDO770:          ENDDO
782:       ENDIF771:       ENDIF
783:       IF (HARMONICF) THEN772:       IF (HARMONICF) THEN
784:          WRITE(MYUNIT,'(A,F12.4)') 'io1> harmonically constrained atoms: strength = ', HARMONICSTR773:          WRITE(MYUNIT,'(A,F12.4)') 'io1> harmonically constrained atoms: strength = ', HARMONICSTR
913: !            WRITE(MYUNIT, '(A,G15.5,A,G15.5)') 'Maximum step size scaled by estimated nearest neighbour distance of ',902: !            WRITE(MYUNIT, '(A,G15.5,A,G15.5)') 'Maximum step size scaled by estimated nearest neighbour distance of ',
914: !    &                    0.677441D0-0.0037582*NATOMS+9.40318D-6*NATOMS**2-6.21931D-9*NATOMS**3,' to give ',STEP(1)903: !    &                    0.677441D0-0.0037582*NATOMS+9.40318D-6*NATOMS**2-6.21931D-9*NATOMS**3,' to give ',STEP(1)
915:          ELSEIF (MULLERBROWNT) THEN 904:          ELSEIF (MULLERBROWNT) THEN 
916:             RADIUS=100.0D0905:             RADIUS=100.0D0
917:          ELSE 906:          ELSE 
918:             RADIUS=RADIUS*2.0D0**(1.0D0/6.0D0)907:             RADIUS=RADIUS*2.0D0**(1.0D0/6.0D0)
919:          ENDIF908:          ENDIF
920:       ENDIF909:       ENDIF
921:       IF ((.NOT.PERIODIC).AND.(.NOT.AMBER).AND.(.NOT.BLNT).AND.(.NOT.MULLERBROWNT).AND.(.NOT.MODEL1T).AND.(.NOT.PERCOLATET) 910:       IF ((.NOT.PERIODIC).AND.(.NOT.AMBER).AND.(.NOT.BLNT).AND.(.NOT.MULLERBROWNT).AND.(.NOT.MODEL1T).AND.(.NOT.PERCOLATET) 
922:      &                    .AND.(.NOT.QCIPOTT).AND.(.NOT.INTCONSTRAINTT).AND.(.NOT.MLP3T).AND.(.NOT.MKTRAPT).AND.(.NOT.MLQT)911:      &                    .AND.(.NOT.QCIPOTT).AND.(.NOT.INTCONSTRAINTT).AND.(.NOT.MLP3T).AND.(.NOT.MKTRAPT).AND.(.NOT.MLQT)
923:      &                    .AND.(.NOT.MLPVB3T).AND.(.NOT.ORBITALS)) 912:      &                    .AND.(.NOT.MLPVB3T)) 
924:      1                    WRITE(MYUNIT,'(A,F20.10)') 'Container radius=',RADIUS913:      1                    WRITE(MYUNIT,'(A,F20.10)') 'Container radius=',RADIUS
925:       RADIUS=RADIUS**2914:       RADIUS=RADIUS**2
926:       IF (PERCOLATET) WRITE(MYUNIT,'(A,F20.10)') 'Checking for percolated structure, cutoff=',PERCCUT915:       IF (PERCOLATET) WRITE(MYUNIT,'(A,F20.10)') 'Checking for percolated structure, cutoff=',PERCCUT
927:       PERCCUT=PERCCUT**2916:       PERCCUT=PERCCUT**2
928:       IF (NPAR.GT.1) THEN917:       IF (NPAR.GT.1) THEN
929:          WRITE(MYUNIT,'(I2,A)') NPAR,' parallel runs'918:          WRITE(MYUNIT,'(I2,A)') NPAR,' parallel runs'
930:          IF (TABOOT) WRITE(MYUNIT,'(A,I4,A)') 'Taboo lists contain the lowest ',NTAB,' minima'919:          IF (TABOOT) WRITE(MYUNIT,'(A,I4,A)') 'Taboo lists contain the lowest ',NTAB,' minima'
931:       ELSE IF (TABOOT) THEN920:       ELSE IF (TABOOT) THEN
932:          WRITE(MYUNIT,'(A,I4,A)') 'Taboo list contains the lowest ',NTAB,' minima'921:          WRITE(MYUNIT,'(A,I4,A)') 'Taboo list contains the lowest ',NTAB,' minima'
933:       ENDIF922:       ENDIF


r32722/keywords.f 2017-06-05 15:30:20.635798367 +0100 r32721/keywords.f 2017-06-05 15:30:22.723825308 +0100
 50: !     &                                 AMBER12_RESIDUES, 50: !     &                                 AMBER12_RESIDUES,
 51: !     &                                 POPULATE_ATOM_DATA 51: !     &                                 POPULATE_ATOM_DATA
 52:       USE CHIRALITY, ONLY : CIS_TRANS_TOL 52:       USE CHIRALITY, ONLY : CIS_TRANS_TOL
 53:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR 53:       USE ISO_C_BINDING, ONLY: C_NULL_CHAR
 54:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS, 54:       USE PARSE_POT_PARAMS, ONLY : PARSE_MGUPTA_PARAMS, PARSE_MSC_PARAMS,
 55:      &     PARSE_MLJ_PARAMS 55:      &     PARSE_MLJ_PARAMS
 56:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT 56:       USE ROTAMER, ONLY: ROTAMER_MOVET, ROTAMER_SCRIPT, ROTAMER_INIT
 57:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE 57:       USE HINGE_MOVES, ONLY: HINGE_INITIALISE
 58:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS 58:       USE MOLECULAR_DYNAMICS, ONLY : MDT, MD_TSTEP, MD_GAMMA, MD_NWAIT, MD_NFREQ, MD_NSTEPS
 59:       USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT 59:       USE OPEP_INTERFACE_MOD, ONLY : OPEP_INIT
 60:       USE ORBITALS_MOD, ONLY: ORBITALS_INIT 
 61:       USE EWALD 60:       USE EWALD
 62:        61:       
 63:       IMPLICIT NONE 62:       IMPLICIT NONE
 64:  63: 
 65:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:) 64:       DOUBLE PRECISION, ALLOCATABLE :: MLPMEAN(:), MLQMEAN(:)
 66:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, J2, J3 65:       INTEGER ITEM, NITEMS, LOC, LINE, NCR, NERROR, LAST, IX, J1, JP, NPCOUNT, NDUMMY, J2, J3
 67:       INTEGER DATA_UNIT, FUNIT 66:       INTEGER DATA_UNIT, FUNIT
 68:       INTEGER MOVABLEATOMINDEX 67:       INTEGER MOVABLEATOMINDEX
 69:       LOGICAL CAT, YESNO, PERMFILE, CONFILE 68:       LOGICAL CAT, YESNO, PERMFILE, CONFILE
 70:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR, 69:       COMMON /BUFINF/ ITEM, NITEMS, LOC(80), LINE, SKIPBL, CLEAR, NCR,
7876:          POLIRT=.TRUE.7875:          POLIRT=.TRUE.
7877: 7876: 
7878:       ELSE IF (WORD.EQ.'MBPOL') THEN7877:       ELSE IF (WORD.EQ.'MBPOL') THEN
7879:          CALL MBPOLINIT7878:          CALL MBPOLINIT
7880:          MBPOLT=.TRUE.7879:          MBPOLT=.TRUE.
7881:          MOLECULART=.TRUE.7880:          MOLECULART=.TRUE.
7882:          ALLOCATE(TYPECH(NATOMSALLOC))7881:          ALLOCATE(TYPECH(NATOMSALLOC))
7883:          TYPECH(1::3)(1:1)="O"7882:          TYPECH(1::3)(1:1)="O"
7884:          TYPECH(2::3)(1:1)="H"7883:          TYPECH(2::3)(1:1)="H"
7885:          TYPECH(3::3)(1:1)="H"7884:          TYPECH(3::3)(1:1)="H"
7886: ! cs675> 
7887:       ELSE IF (WORD.EQ.'ORBITALS') THEN 
7888:          ORBITALS=.TRUE. 
7889:          CALL READI(ORBVAREXPONENT) 
7890:          CALL ORBITALS_INIT(NORBS, NROTS) 
7891:       ELSE7885:       ELSE
7892:          CALL REPORT('Unrecognized command '//WORD,.TRUE.)7886:          CALL REPORT('Unrecognized command '//WORD,.TRUE.)
7893:          STOP7887:          STOP
7894:       ENDIF7888:       ENDIF
7895:       CALL FLUSH(MYUNIT)7889:       CALL FLUSH(MYUNIT)
7896: 7890: 
7897: !ds656> NTYPEA can fluctuate in homotop refinemet routine, so we need7891: !ds656> NTYPEA can fluctuate in homotop refinemet routine, so we need
7898: !     a fixed reference. It is set here instead of every single7892: !     a fixed reference. It is set here instead of every single
7899: !     IF block where NTYPEA can potentially change.7893: !     IF block where NTYPEA can potentially change.
7900:       NTYPEA_FIX = NTYPEA7894:       NTYPEA_FIX = NTYPEA


r32722/main.F 2017-06-05 15:30:20.863801306 +0100 r32721/main.F 2017-06-05 15:30:22.943828146 +0100
131:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS,131:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS,
132:      &                         ' should be ',MLPHIDDEN*(MLPIN+MLPOUT) +1132:      &                         ' should be ',MLPHIDDEN*(MLPIN+MLPOUT) +1
133:             STOP133:             STOP
134:          ENDIF134:          ENDIF
135:       ELSEIF (MLP3T) THEN135:       ELSEIF (MLP3T) THEN
136:          IF (NUMBER_OF_ATOMS.NE.MLPHIDDEN*(MLPIN+MLPOUT)) THEN136:          IF (NUMBER_OF_ATOMS.NE.MLPHIDDEN*(MLPIN+MLPOUT)) THEN
137:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS,137:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS,
138:      &                         ' should be ',MLPHIDDEN*(MLPIN+MLPOUT) 138:      &                         ' should be ',MLPHIDDEN*(MLPIN+MLPOUT) 
139:             STOP139:             STOP
140:          ENDIF140:          ENDIF
141:       ELSEIF (ORBITALS) THEN 
142:          IF (NUMBER_OF_ATOMS.NE.NROTS) THEN 
143:             WRITE(MYUNIT,'(A,I8,A,I8)') 'main> ERROR *** number of coordinates in coords is ',NUMBER_OF_ATOMS, 
144:      &                         ' should be ',NROTS 
145:             STOP 
146:          ENDIF 
147:       ENDIF141:       ENDIF
148:       NATOMSALLOC=NUMBER_OF_ATOMS142:       NATOMSALLOC=NUMBER_OF_ATOMS
149:       NATOMS=NUMBER_OF_ATOMS143:       NATOMS=NUMBER_OF_ATOMS
150:       IF (GCBHT) THEN144:       IF (GCBHT) THEN
151:          NATOMSALLOC=GCNATOMS145:          NATOMSALLOC=GCNATOMS
152:          IF (GCNATOMS.LT.NUMBER_OF_ATOMS) THEN146:          IF (GCNATOMS.LT.NUMBER_OF_ATOMS) THEN
153:             WRITE(MYUNIT,'(A)') 'main> *** ERROR - number of atoms exceeds maximum allowed by GCBH keyword'147:             WRITE(MYUNIT,'(A)') 'main> *** ERROR - number of atoms exceeds maximum allowed by GCBH keyword'
154:             STOP148:             STOP
155:          ENDIF149:          ENDIF
156:          ALLOCATE(QENERGIES(GCNATOMS),QCOORDINATES(GCNATOMS,3*GCNATOMS),QPE(GCNATOMS))150:          ALLOCATE(QENERGIES(GCNATOMS),QCOORDINATES(GCNATOMS,3*GCNATOMS),QPE(GCNATOMS))


r32722/orbitals.f90 2017-06-05 15:30:21.079804094 +0100 r32721/orbitals.f90 2017-06-05 15:30:23.163830986 +0100
  1: MODULE ORBITALS_MOD  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/GMIN/source/orbitals.f90' in revision 32721
  2: IMPLICIT NONE 
  3:  
  4:   CONTAINS 
  5:  
  6:     SUBROUTINE ORBITALS_INIT(NORBS, NROTS) 
  7:  
  8:         ! Initialise information needed for orbital optimisation. 
  9:         ! This requires reading in the integral arrays from file. 
 10:  
 11:         USE COMMONS, ONLY: R2INTS, DIPINTS, ORBVAREXPONENT 
 12:  
 13:         INTEGER, INTENT(OUT) :: NORBS, NROTS 
 14:  
 15:         INTEGER :: ORB_FILE_UNIT 
 16:  
 17:         CHARACTER (LEN=10) :: WORD 
 18:  
 19:         LOGICAL :: END 
 20:  
 21:         CALL FILE_OPEN('INTEGRALS_INFO',ORB_FILE_UNIT,.FALSE.) 
 22:         CALL INPUT(END, ORB_FILE_UNIT) 
 23:  
 24:         CALL READI(NORBS) 
 25:         NROTS = NORBS * (NORBS - 1) / 2 
 26:  
 27:         !PRINT '(A,3I10)','NORBS,NROTS,ORBVAREXPONENT=',NORBS,NROTS,ORBVAREXPONENT 
 28:  
 29:         CLOSE(ORB_FILE_UNIT) 
 30:  
 31:         ALLOCATE(R2INTS(NORBS,NORBS)) 
 32:         ALLOCATE(DIPINTS(3,NORBS,NORBS)) 
 33:  
 34:         !PRINT '(A)','R2' 
 35:         CALL READ_INTEGRALS('INTEGRALS_R2.csv', NORBS, R2INTS(:,:)) 
 36:         !PRINT '(A)','X' 
 37:         CALL READ_INTEGRALS('INTEGRALS_X.csv', NORBS, DIPINTS(1,:,:)) 
 38:         !PRINT '(A)','Y' 
 39:         CALL READ_INTEGRALS('INTEGRALS_Y.csv', NORBS, DIPINTS(2,:,:)) 
 40:         !PRINT '(A)','Z' 
 41:         CALL READ_INTEGRALS('INTEGRALS_Z.csv', NORBS, DIPINTS(3,:,:)) 
 42:  
 43:  
 44:     END SUBROUTINE ORBITALS_INIT 
 45:  
 46:     SUBROUTINE READ_INTEGRALS(INT_FILENAME, NORBS, INTS) 
 47:  
 48:         ! Given filename of intergral files reads the integrals within 
 49:         ! it into the array INTS. Integrals are assumed to always have 
 50:         ! dimension NORBSxNORBS. 
 51:  
 52:         CHARACTER (LEN=*), INTENT(IN) :: INT_FILENAME 
 53:         INTEGER, INTENT(IN) :: NORBS 
 54:         DOUBLE PRECISION, INTENT(OUT) :: INTS(NORBS, NORBS) 
 55:         CHARACTER (LEN=32) :: CALCID 
 56:         INTEGER :: INT_FILE_UNIT, I 
 57:         LOGICAL :: END 
 58:  
 59:         CALL FILE_OPEN(INT_FILENAME,INT_FILE_UNIT,.FALSE.) 
 60:         CALL INPUT(END, INT_FILE_UNIT) 
 61:         CALL READU(CALCID) 
 62:         !PRINT '(A,I10,X,A)','INT_FILE_UNIT,CALCID=',INT_FILE_UNIT,CALCID 
 63:         !PRINT '(A,I10,A)','NORBS=',NORBS 
 64:  
 65:         DO I = 1, NORBS 
 66:            READ(INT_FILE_UNIT,*) INTS(I,1:NORBS)  
 67:            !PRINT *,INTS(I,1:NORBS) 
 68:         END DO 
 69:  
 70:     END SUBROUTINE READ_INTEGRALS 
 71:  
 72:     SUBROUTINE GET_ORBITAL_LOCALITY(X, GRAD, LOCALITY, GTEST) 
 73:  
 74:         ! Obtains the Power of the Orbital Variance locality measure for 
 75:         ! a given position. If GTEST is set also returns the gradient. 
 76:  
 77:         USE COMMONS, ONLY: NROTS, NORBS, R2INTS, DIPINTS, ORBVAREXPONENT 
 78:  
 79:         DOUBLE PRECISION, INTENT(IN) :: X(NROTS) 
 80:         DOUBLE PRECISION, INTENT(OUT) :: GRAD(NROTS), LOCALITY 
 81:         LOGICAL, INTENT(IN) :: GTEST 
 82:  
 83:         DOUBLE PRECISION :: ROTATION(NORBS,NORBS), ORBVAR(NORBS) 
 84:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS), IND_ORBVAR_DERIV(NORBS,NORBS) 
 85:         DOUBLE PRECISION :: PENALTY_DERIV(NORBS,NORBS), ROT_DERIV(NORBS,NORBS) 
 86:         DOUBLE PRECISION :: PROD_DERIVS(NORBS,NORBS) 
 87:  
 88:         INTEGER :: I,J,K 
 89:  
 90:         CALL GET_ROTATION(X, ROTATION) 
 91:  
 92:         CALL GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES) 
 93:         DO I = 1, NORBS 
 94:             ORBVAR(I) = -SUM(ORBDIPOLES(1:3,I)**2) 
 95:             DO J = 1, NORBS 
 96:                 DO K = 1, NORBS 
 97:                     ORBVAR(I)=ORBVAR(I)+R2INTS(J,K)*ROTATION(I,J)*ROTATION(I,K) 
 98:                 END DO 
 99:             END DO 
100:             !PRINT *,'ORBITAL',I,',ORBVAR=',ORBVAR(I),',DIPOLES=',ORBDIPOLES(1:3,I) 
101:  
102:         END DO 
103:  
104:  
105:         !PRINT *,'ORBVARS=',ORBVAR 
106:         !PRINT *,'ORBVARSEXP=',ORBVAR**ORBVAREXPONENT 
107:  
108:         LOCALITY = SUM(ORBVAR ** ORBVAREXPONENT) 
109:  
110:         !PRINT *,'LOCALITY=',LOCALITY 
111:  
112:         IF (GTEST) THEN 
113:             !PRINT *,'ORBDIPOLES=' 
114:             !DO J = 1,3 
115:             !    PRINT *,ORBDIPOLES(J,1:NORBS) 
116:             !END DO 
117:  
118:             CALL GET_INDIVIDUAL_DERIVATIVES(ROTATION, ORBDIPOLES, IND_ORBVAR_DERIV) 
119:             !PRINT *,'IND_ORBVAR_DERIV=' 
120:             !DO J = 1,NORBS 
121:             !    PRINT *, IND_ORBVAR_DERIV(J,1:NORBS) 
122:             !END DO 
123:  
124:  
125:             CALL GET_PENALTY_DERIVATIVES(IND_ORBVAR_DERIV, ORBVAR, PENALTY_DERIV) 
126:  
127:             !PRINT *,'PEN_DERIV=' 
128:             !DO J = 1,NORBS 
129:             !    PRINT *, PENALTY_DERIV(J,1:NORBS) 
130:             !END DO 
131:  
132:             DO I = 1, NROTS 
133:                 CALL GET_ROTATION_DERIVATIVE(X, I, ROT_DERIV) 
134:                 !PRINT *,'ROT_DERIV=' 
135:                 !DO J = 1,NORBS 
136:                 !    PRINT *, ROT_DERIV(J,1:NORBS) 
137:                 !END DO 
138:  
139:                 CALL ELEMENTWISE_MULTIPLICATION(ROT_DERIV,PENALTY_DERIV,PROD_DERIVS) 
140:                 !PRINT *,'OVERALL_PROD=' 
141:                 !DO J = 1,NORBS 
142:                 !    PRINT *, PROD_DERIVS(J,:) 
143:                 !END DO 
144:                 GRAD(I) = SUM(PROD_DERIVS) 
145:                 !PRINT *,'GRAD',I,'=',GRAD(I) 
146:             END DO 
147:  
148:         END IF 
149:  
150:     END SUBROUTINE GET_ORBITAL_LOCALITY 
151:  
152:     SUBROUTINE ELEMENTWISE_MULTIPLICATION(MATA,MATB,MATC) 
153:  
154:         USE COMMONS, ONLY: NORBS 
155:  
156:         DOUBLE PRECISION, INTENT(IN) :: MATA(NORBS,NORBS), MATB(NORBS,NORBS) 
157:  
158:         DOUBLE PRECISION, INTENT(OUT) :: MATC(NORBS, NORBS) 
159:  
160:         INTEGER :: I, J 
161:  
162:         DO I = 1, NORBS 
163:             DO J = 1, NORBS 
164:                 MATC(I,J) = MATA(I,J) * MATB(I,J) 
165:             END DO 
166:         END DO 
167:  
168:     END SUBROUTINE ELEMENTWISE_MULTIPLICATION 
169:  
170:     SUBROUTINE GET_ROTATION(COORDS, ROTATION) 
171:  
172:         ! Obtains the overall rotation resulting from the direct product 
173:         ! of the Givens rotations corresponding to the specified coordinates. 
174:  
175:         USE COMMONS, ONLY: NROTS, NORBS 
176:  
177:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
178:  
179:         DOUBLE PRECISION, INTENT(OUT) :: ROTATION(NORBS,NORBS) 
180:         DOUBLE PRECISION :: NEXT_ROT(NORBS,NORBS) 
181:  
182:         INTEGER :: I, J, TRIIND 
183:  
184:         ROTATION(:,:) = 0.0 
185:  
186:         DO TRIIND = 1, NORBS 
187:             ROTATION(TRIIND, TRIIND) = 1.0 
188:         END DO 
189:  
190:         DO TRIIND = 1, NROTS 
191:             CALL DECODE_TRIIND(TRIIND, I, J) 
192:             CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
193:             ROTATION = MATMUL(ROTATION,NEXT_ROT) 
194:         END DO 
195:  
196:     END SUBROUTINE GET_ROTATION 
197:  
198:     SUBROUTINE GET_ROTATION_DERIVATIVE(COORDS, DERIVTRIIND, DERIV_ROTATION) 
199:  
200:         ! Obtain the derivative of the rotation matrix with respect to a single 
201:         ! rotation angle theta at the current coordinates. 
202:  
203:         USE COMMONS, ONLY: NROTS, NORBS 
204:  
205:         DOUBLE PRECISION, INTENT(IN) :: COORDS(NROTS) 
206:         INTEGER, INTENT(IN) :: DERIVTRIIND 
207:  
208:         DOUBLE PRECISION, INTENT(OUT) :: DERIV_ROTATION(NORBS,NORBS) 
209:  
210:         DOUBLE PRECISION :: NEXT_ROT(NORBS,NORBS) 
211:         INTEGER :: TRIIND, I, J 
212:  
213:         DERIV_ROTATION(:,:) = 0.0 
214:  
215:         DO TRIIND = 1, NORBS 
216:             DERIV_ROTATION(TRIIND, TRIIND) = 1.0 
217:         END DO 
218:  
219:         DO TRIIND = 1, NROTS 
220:             CALL DECODE_TRIIND(TRIIND, I, J) 
221:             IF (TRIIND.EQ.DERIVTRIIND) THEN 
222:                 CALL GET_DERIVATIVE_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
223:             ELSE 
224:                 CALL GET_GIVENS_ROTATION(I,J,COORDS(TRIIND),NEXT_ROT) 
225:             END IF 
226:             DERIV_ROTATION = MATMUL(DERIV_ROTATION,NEXT_ROT) 
227:         END DO 
228:  
229:         ! This appears to be required- investigate why. Something confusing in 
230:         ! fortran array indexing. 
231:         DERIV_ROTATION = TRANSPOSE(DERIV_ROTATION) 
232:         !DO TRIIND = 1, NORBS 
233:         !    PRINT *,TRIIND,'DERIV_ROT=',DERIV_ROTATION(TRIIND,1:NORBS) 
234:         !END DO 
235:  
236:     END SUBROUTINE GET_ROTATION_DERIVATIVE 
237:  
238:     SUBROUTINE GET_GIVENS_ROTATION(I, J, THETA, MAT) 
239:  
240:         ! Obtains Givens rotation between orbitals I and J with angle theta. 
241:  
242:         USE COMMONS, ONLY: NORBS 
243:  
244:         INTEGER, INTENT(IN) :: I, J 
245:         DOUBLE PRECISION, INTENT(IN) :: THETA 
246:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS) 
247:  
248:         INTEGER :: VAL 
249:  
250:         MAT(1:NORBS,1:NORBS) = 0.0 
251:  
252:         DO VAL = 1, NORBS 
253:             MAT(VAL, VAL) = 1.0 
254:         END DO 
255:  
256:         MAT(I,I) = COS(THETA) 
257:         MAT(J,J) = COS(THETA) 
258:         MAT(I,J) = -SIN(THETA) 
259:         MAT(J,I) = SIN(THETA) 
260:  
261:         !PRINT *, 'GIVENSROT',I,'-',J,'=' 
262:         !DO VAL=1,NORBS 
263:         !    PRINT *,MAT(VAL,1:NORBS) 
264:         !END DO 
265:  
266:     END SUBROUTINE GET_GIVENS_ROTATION 
267:  
268:     SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION(I, J, THETA, MAT) 
269:  
270:         ! Obtains the derivative of the Givens rotation between orbitals I and J with angle 
271:         ! theta as a matrix wrt the rotation angle theta. 
272:  
273:         USE COMMONS, ONLY: NORBS 
274:  
275:         INTEGER, INTENT(IN) :: I, J 
276:         DOUBLE PRECISION, INTENT(IN) :: THETA 
277:         DOUBLE PRECISION, INTENT(OUT) :: MAT(NORBS,NORBS) 
278:  
279:         INTEGER :: VAL 
280:  
281:         MAT(1:NORBS,1:NORBS) = 0.0 
282:  
283:         MAT(I,I) = -SIN(THETA) 
284:         MAT(J,J) = -SIN(THETA) 
285:         MAT(I,J) = -COS(THETA) 
286:         MAT(J,I) = COS(THETA) 
287:         !PRINT *, 'DERIVGIVENSROT',I,'-',J,'=' 
288:         !DO VAL=1,NORBS 
289:         !    PRINT *,MAT(VAL,1:NORBS) 
290:         !END DO 
291:  
292:     END SUBROUTINE GET_DERIVATIVE_GIVENS_ROTATION 
293:  
294:     SUBROUTINE GET_CURRENT_ORB_DIPOLES(ROTATION, ORBDIPOLES) 
295:  
296:         ! Obtain the dipole moments of the orbital set created by the current rotation of the starting MOs. 
297:  
298:         USE COMMONS, ONLY: NORBS, DIPINTS 
299:  
300:         DOUBLE PRECISION, INTENT(IN) :: ROTATION(NORBS, NORBS) 
301:         DOUBLE PRECISION :: ORBDIPOLES(3,NORBS) 
302:  
303:         INTEGER :: CURRENT_ORB, ORIG_ORBA, ORIG_ORBB, DIRECTION 
304:  
305:         ORBDIPOLES(1:3,1:NORBS) = 0.0 
306:  
307:         DO DIRECTION = 1,3 
308:             DO CURRENT_ORB = 1, NORBS 
309:                 DO ORIG_ORBA = 1, NORBS 
310:                     DO ORIG_ORBB = 1, NORBS 
311:                         ORBDIPOLES(DIRECTION,CURRENT_ORB) = ORBDIPOLES(DIRECTION,CURRENT_ORB) + & 
312:                             ROTATION(CURRENT_ORB,ORIG_ORBA) * ROTATION(CURRENT_ORB,ORIG_ORBB) * & 
313:                             DIPINTS(DIRECTION,ORIG_ORBA,ORIG_ORBB) 
314:                     END DO 
315:                 END DO 
316:             END DO 
317:         END DO 
318:  
319:     END SUBROUTINE GET_CURRENT_ORB_DIPOLES 
320:  
321:     SUBROUTINE DECODE_TRIIND(TRIIND, I, J) 
322:  
323:         ! Given a triangular index obtains the 2D indices corresponding to it. 
324:  
325:         USE COMMONS, ONLY: NORBS 
326:  
327:         INTEGER, INTENT(IN) :: TRIIND 
328:         INTEGER, INTENT(OUT) :: I, J 
329:  
330:         INTEGER :: RUNNING_SUM, VAL, I_, J_ 
331:  
332:         RUNNING_SUM = 0 
333:  
334:         DO VAL = 1, NORBS 
335:             IF (RUNNING_SUM + VAL .GE. TRIIND) THEN 
336:                 I_ = VAL + 1 
337:                 J_ = TRIIND - RUNNING_SUM 
338:                 EXIT 
339:             ELSE 
340:                 RUNNING_SUM = RUNNING_SUM + VAL 
341:             END IF 
342:  
343:         END DO 
344:  
345:         IF (I_ < J_) THEN 
346:             I = J_ 
347:             J = I_ 
348:         ELSE 
349:             I = I_ 
350:             J = J_ 
351:         END IF 
352:  
353:         !PRINT '(A,3I10)','TRIIND,I,J=',TRIIND,I,J 
354:  
355:     END SUBROUTINE DECODE_TRIIND 
356:  
357:     SUBROUTINE GET_INDIVIDUAL_DERIVATIVES(CURRENT_ROTATION, CURRENT_DIP, INDIVIDUAL_DERIVATIVES) 
358:  
359:         ! Generates the matrix of all derivatives of different orbital variances w.r.t. 
360:         ! the coefficients of their expansion in the original MO basis. 
361:  
362:         USE COMMONS, ONLY: NORBS 
363:  
364:         DOUBLE PRECISION, INTENT(IN) :: CURRENT_ROTATION(NORBS, NORBS), CURRENT_DIP(3,NORBS) 
365:  
366:         DOUBLE PRECISION, INTENT(OUT) :: INDIVIDUAL_DERIVATIVES(NORBS,NORBS) 
367:  
368:         INTEGER :: ORIGORB, NEWORB 
369:  
370:         !PRINT *,'CURRENTROTATION=',CURRENT_ROTATION 
371:  
372:         DO NEWORB = 1, NORBS 
373:             DO ORIGORB = 1, NORBS 
374:                 CALL GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEWORB, ORIGORB, INDIVIDUAL_DERIVATIVES(NEWORB,ORIGORB)) 
375:             END DO 
376:         END DO 
377:  
378:     END SUBROUTINE GET_INDIVIDUAL_DERIVATIVES 
379:  
380:     SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE(CURRENT_ROTATION, CURRENT_DIP, NEW_ORB, ORIG_ORB, DERIVATIVE) 
381:  
382:         ! Obtains the derivative of a single orbital at the current rotation's orbital 
383:         ! variance w.r.t. to the coefficient an original (unrotated) orbital. 
384:  
385:         USE COMMONS, ONLY: NORBS, R2INTS, DIPINTS 
386:  
387:         DOUBLE PRECISION, INTENT(IN) :: CURRENT_ROTATION(NORBS, NORBS), CURRENT_DIP(3,NORBS) 
388:  
389:         DOUBLE PRECISION, INTENT(OUT) :: DERIVATIVE 
390:  
391:         INTEGER,INTENT(IN) :: NEW_ORB, ORIG_ORB 
392:  
393:         INTEGER :: OTHER_ORIG_ORB 
394:  
395:         DOUBLE PRECISION :: CONTRIB 
396:  
397:         DERIVATIVE = 0 
398:  
399:         DO OTHER_ORIG_ORB = 1, NORBS 
400:             CONTRIB = 2 * CURRENT_ROTATION(NEW_ORB,OTHER_ORIG_ORB)*R2INTS(ORIG_ORB, OTHER_ORIG_ORB) 
401:  
402:             DERIVATIVE = DERIVATIVE + CONTRIB 
403:  
404:  
405:             CONTRIB = 4 * CURRENT_ROTATION(NEW_ORB,OTHER_ORIG_ORB)*& 
406:                         (CURRENT_DIP(1,NEW_ORB)*DIPINTS(1,ORIG_ORB,OTHER_ORIG_ORB)& 
407:                         +CURRENT_DIP(2,NEW_ORB)*DIPINTS(2,ORIG_ORB,OTHER_ORIG_ORB)& 
408:                         +CURRENT_DIP(3,NEW_ORB)*DIPINTS(3,ORIG_ORB,OTHER_ORIG_ORB)) 
409:  
410:             DERIVATIVE = DERIVATIVE - CONTRIB 
411:  
412:         END DO 
413:  
414:     END SUBROUTINE GET_SINGLE_ORBVAR_DERIVATIVE 
415:  
416:     SUBROUTINE GET_PENALTY_DERIVATIVES(IND_ORBVAR_DERIV, ORB_VAR, PENALTY_DERIVATIVES) 
417:  
418:         ! Get the derivative of the overall penalty function with respect to each element 
419:         ! of the rotation matrix. 
420:  
421:         USE COMMONS, ONLY: NORBS, ORBVAREXPONENT 
422:  
423:         DOUBLE PRECISION, INTENT(IN) :: IND_ORBVAR_DERIV(NORBS,NORBS), ORB_VAR(NORBS) 
424:         DOUBLE PRECISION, INTENT(OUT) :: PENALTY_DERIVATIVES(NORBS, NORBS) 
425:  
426:         INTEGER :: ORIG_ORB, NEW_ORB 
427:  
428:         PENALTY_DERIVATIVES = 0.0 
429:  
430:         DO ORIG_ORB = 1, NORBS 
431:             DO NEW_ORB = 1, NORBS 
432:                 PENALTY_DERIVATIVES(ORIG_ORB, NEW_ORB) = & 
433:                     ORBVAREXPONENT * IND_ORBVAR_DERIV(NEW_ORB,ORIG_ORB) * ORB_VAR(NEW_ORB) ** (ORBVAREXPONENT-1) 
434:             END DO 
435:         END DO 
436:  
437:     END SUBROUTINE GET_PENALTY_DERIVATIVES 
438:  
439:  
440:     !SUBROUTINE ORBITALS_FINISH() 
441:  
442:         ! Deallocate arrays in commons? 
443:  
444:     !END SUBROUTINE ORBITALS_FINISH 
445: END MODULE ORBITALS_MOD 


r32722/potential.f90 2017-06-05 15:30:21.299806932 +0100 r32721/potential.f90 2017-06-05 15:30:23.387833876 +0100
 15: !   You should have received a copy of the GNU General Public License 15: !   You should have received a copy of the GNU General Public License
 16: !   along with this program; if not, write to the Free Software 16: !   along with this program; if not, write to the Free Software
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19: SUBROUTINE POTENTIAL(X, GRAD, EREAL, GRADT, SECT) 19: SUBROUTINE POTENTIAL(X, GRAD, EREAL, GRADT, SECT)
 20: ! This subroutine acts as a wrapper for all of the potentials implemented in GMIN. 20: ! This subroutine acts as a wrapper for all of the potentials implemented in GMIN.
 21: ! The potential is chosen in the 'data' file (default is Lennard-Jones particles). 21: ! The potential is chosen in the 'data' file (default is Lennard-Jones particles).
 22: ! 22: !
 23: ! In addition, this routine also applies fields, if appropriate. 23: ! In addition, this routine also applies fields, if appropriate.
 24:    USE COMMONS 24:    USE COMMONS
 25:    USE ORBITALS_MOD 
 26:    USE GENRIGID, ONLY: RIGIDINIT, NRIGIDBODY, DEGFREEDOMS, ATOMRIGIDCOORDT, & 25:    USE GENRIGID, ONLY: RIGIDINIT, NRIGIDBODY, DEGFREEDOMS, ATOMRIGIDCOORDT, &
 27:                        AACONVERGENCET, AACONVERGENCE, TRANSFORMRIGIDTOC, & 26:                        AACONVERGENCET, AACONVERGENCE, TRANSFORMRIGIDTOC, &
 28:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, & 27:                        TRANSFORMGRAD, TRANSFORMHESSIAN, TRANSFORMCTORIGID, &
 29:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS, aaconvergence2 28:                        GENRIGID_DISTANCECHECK, GENRIGID_COMPRESS, aaconvergence2
 30:    USE MULTIPOT, ONLY: MULTIPOT_CALL 29:    USE MULTIPOT, ONLY: MULTIPOT_CALL
 31:    USE QMODULE, ONLY: QMIN, QMINP 30:    USE QMODULE, ONLY: QMIN, QMINP
 32:    USE PERMU, ONLY: FIN 31:    USE PERMU, ONLY: FIN
 33:    USE PORFUNCS 32:    USE PORFUNCS
 34:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C,& 33:    USE AMBER12_INTERFACE_MOD, ONLY: AMBER12_ENERGY_AND_GRADIENT, POT_ENE_REC_C,&
 35: &                                   AMBER12_NUM_HESS 34: &                                   AMBER12_NUM_HESS
421: !          X(J1)=X(J1)+DIFF420: !          X(J1)=X(J1)+DIFF
422: ! !        IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN421: ! !        IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN
423: !             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)422: !             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF)
424: ! !        ENDIF423: ! !        ENDIF
425: !       ENDDO424: !       ENDDO
426: 425: 
427:    ELSE IF (MLP3T) THEN426:    ELSE IF (MLP3T) THEN
428:       CALL MLP3(X, GRAD, EREAL, GRADT, SECT)427:       CALL MLP3(X, GRAD, EREAL, GRADT, SECT)
429:    ELSEIF (MLQT) THEN428:    ELSEIF (MLQT) THEN
430:       CALL MLQ(X, GRAD, EREAL, GRADT, SECT)429:       CALL MLQ(X, GRAD, EREAL, GRADT, SECT)
431:    ELSE IF (ORBITALS) THEN 
432:       CALL GET_ORBITAL_LOCALITY(X, GRAD, EREAL, GRADT) 
433:       DIFF=1.0D-4 
434:       WRITE(MYUNIT, *) 'analytic and numerical gradients:' 
435:       DO J1=1, NATOMS 
436:          X(J1)=X(J1)+DIFF 
437:          CALL GET_ORBITAL_LOCALITY(X, GPLUS, EPLUS,.FALSE.) 
438:          X(J1)=X(J1)-2.0D0*DIFF 
439:          CALL GET_ORBITAL_LOCALITY(X, GMINUS, EMINUS,.FALSE.) 
440:          X(J1)=X(J1)+DIFF 
441:          IF ((ABS(GRAD(J1)).NE.0.0D0).AND.(100.0D0*(GRAD(J1)-(EPLUS-EMINUS)/(2.0D0*DIFF))/GRAD(J1).GT.1.0D0)) THEN 
442:             WRITE(MYUNIT,'(A,I5, 2F20.10)') 'gtest ', J1, GRAD(J1),(EPLUS-EMINUS)/(2.0D0*DIFF) 
443:          ENDIF 
444:       ENDDO 
445:       STOP 
446:    ELSE IF (LJATT) THEN430:    ELSE IF (LJATT) THEN
447:       CALL LJ(X, GRAD, EREAL, GRADT, SECT)431:       CALL LJ(X, GRAD, EREAL, GRADT, SECT)
448:       CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)432:       CALL AXT(NATOMS, X, GRAD, EREAL, GRADT, ZSTAR)
449: 433: 
450:    ELSE IF (DFTBCT) THEN434:    ELSE IF (DFTBCT) THEN
451:       IF (.NOT.PERCOLATET) CALL RAD(X, GRAD, EREAL, GRADT)435:       IF (.NOT.PERCOLATET) CALL RAD(X, GRAD, EREAL, GRADT)
452:       CALL DFTBC(NATOMS, X, GRAD, EREAL, GRADT)436:       CALL DFTBC(NATOMS, X, GRAD, EREAL, GRADT)
453:       IF (FTEST) THEN437:       IF (FTEST) THEN
454:          RETURN438:          RETURN
455:       END IF439:       END IF
1529:             CALL HARMONICFIELD( X(3*J1-2), HARMONICR0(3*J1-2), GRAD(3*J1-2), EREAL)1513:             CALL HARMONICFIELD( X(3*J1-2), HARMONICR0(3*J1-2), GRAD(3*J1-2), EREAL)
1530: !            WRITE(*,*) "#V(1) after ", GRAD((J1-1)*3+1)1514: !            WRITE(*,*) "#V(1) after ", GRAD((J1-1)*3+1)
1531:          END IF1515:          END IF
1532:       END DO1516:       END DO
1533:    END IF1517:    END IF
1534: 1518: 
1535:    IF (GRADT .OR. CSMT) THEN1519:    IF (GRADT .OR. CSMT) THEN
1536:       IF (FREEZE) THEN1520:       IF (FREEZE) THEN
1537:          DO J1=1, NATOMS1521:          DO J1=1, NATOMS
1538:             IF (FROZEN(J1)) THEN1522:             IF (FROZEN(J1)) THEN
1539:                IF (MLPB3T.OR.MLQT.OR.MLPVB3T.OR.ORBITALS) THEN1523:                IF (MLPB3T.OR.MLQT.OR.MLPVB3T) THEN
1540:                   GRAD(J1)=0.0D01524:                   GRAD(J1)=0.0D0
1541:                ELSE1525:                ELSE
1542:                   GRAD(3*J1-2 : 3*J1)=0.0D01526:                   GRAD(3*J1-2 : 3*J1)=0.0D0
1543:                ENDIF1527:                ENDIF
1544:             END IF1528:             END IF
1545:          END DO1529:          END DO
1546:       END IF1530:       END IF
1547: 1531: 
1548:       IF (SEEDT .AND. FREEZECORE) THEN1532:       IF (SEEDT .AND. FREEZECORE) THEN
1549:          GRAD(3*(NATOMS-NSEED)+1 : 3*NATOMS)=0.0D01533:          GRAD(3*(NATOMS-NSEED)+1 : 3*NATOMS)=0.0D0
1573: 1557: 
1574:       IF (CSMT .AND. (.NOT.SYMMETRIZECSM)) THEN1558:       IF (CSMT .AND. (.NOT.SYMMETRIZECSM)) THEN
1575:          DUMMY2=0.0D01559:          DUMMY2=0.0D0
1576:          RMS=0.0D01560:          RMS=0.0D0
1577:       ELSE IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN1561:       ELSE IF (RIGIDINIT .AND. (.NOT. ATOMRIGIDCOORDT)) THEN
1578:          DUMMY2=SUM(GRAD(1:DEGFREEDOMS)**2)1562:          DUMMY2=SUM(GRAD(1:DEGFREEDOMS)**2)
1579:          RMS=MAX(SQRT(DUMMY2/DEGFREEDOMS), 1.0D-100)1563:          RMS=MAX(SQRT(DUMMY2/DEGFREEDOMS), 1.0D-100)
1580:          IF ((RMS < 5.0D0 * BQMAX) .AND. AACONVERGENCET .AND. (.NOT. COMPRESSRIGIDT) .AND. (.NOT. COMPON)) THEN1564:          IF ((RMS < 5.0D0 * BQMAX) .AND. AACONVERGENCET .AND. (.NOT. COMPRESSRIGIDT) .AND. (.NOT. COMPON)) THEN
1581:             CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, RMS)1565:             CALL AACONVERGENCE (GRADATOMS, XRIGIDCOORDS, XRIGIDGRAD, RMS)
1582:          END IF1566:          END IF
1583:       ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT.OR.ORBITALS) THEN1567:       ELSE IF (MLP3T.OR.MLPB3T.OR.MLQT) THEN
1584:          DUMMY2=SUM(GRAD(1:NATOMS)**2)1568:          DUMMY2=SUM(GRAD(1:NATOMS)**2)
1585:          RMS=MAX(DSQRT(DUMMY2/(NATOMS)), 1.0D-100)1569:          RMS=MAX(DSQRT(DUMMY2/(NATOMS)), 1.0D-100)
1586:       ELSE IF (.NOT.THOMSONT) THEN1570:       ELSE IF (.NOT.THOMSONT) THEN
1587:          DUMMY2=SUM(GRAD(1:3*NATOMS)**2)1571:          DUMMY2=SUM(GRAD(1:3*NATOMS)**2)
1588:          RMS=MAX(DSQRT(DUMMY2/(3*NATOMS)), 1.0D-100)1572:          RMS=MAX(DSQRT(DUMMY2/(3*NATOMS)), 1.0D-100)
1589:       ELSE1573:       ELSE
1590:          DUMMY2=SUM(GRAD(1:2*NATOMS)**2)1574:          DUMMY2=SUM(GRAD(1:2*NATOMS)**2)
1591:          RMS=MAX(DSQRT(DUMMY2/(2*NATOMS)), 1.0D-100)1575:          RMS=MAX(DSQRT(DUMMY2/(2*NATOMS)), 1.0D-100)
1592:       END IF1576:       END IF
1593:       IF(DEBUG.AND.(RMS.NE.RMS)) THEN1577:       IF(DEBUG.AND.(RMS.NE.RMS)) THEN


r32722/quench.F 2017-06-05 15:30:21.519809772 +0100 r32721/quench.F 2017-06-05 15:30:23.627836971 +0100
109:       IF (DFTBCT.AND.GUIDET) THEN109:       IF (DFTBCT.AND.GUIDET) THEN
110:          LJATT=.TRUE.110:          LJATT=.TRUE.
111:          IF (DEBUG) WRITE(MYUNIT,'(A)') 'quench> Turning on LJAT guiding potential and rescaling coordinates'111:          IF (DEBUG) WRITE(MYUNIT,'(A)') 'quench> Turning on LJAT guiding potential and rescaling coordinates'
112:          COORDS(1:3*NATOMS,NP)=COORDS(1:3*NATOMS,NP)/LJATTOC112:          COORDS(1:3*NATOMS,NP)=COORDS(1:3*NATOMS,NP)/LJATTOC
113:       ENDIF113:       ENDIF
114:       IF (CSMGUIDET) CSMDOGUIDET=.TRUE.114:       IF (CSMGUIDET) CSMDOGUIDET=.TRUE.
115:       NOPT=3*NATOMS115:       NOPT=3*NATOMS
116:       IF (WENZEL) NOPT=2116:       IF (WENZEL) NOPT=2
117:       IF (MULLERBROWNT) NOPT=2117:       IF (MULLERBROWNT) NOPT=2
118:       IF (MLP3T.OR.MLPB3T.OR.MLPVB3T) NOPT=NMLP118:       IF (MLP3T.OR.MLPB3T.OR.MLPVB3T) NOPT=NMLP
119:       IF (ORBITALS) NOPT=NROTS 
120:       IF (MLQT) NOPT=NMLQ119:       IF (MLQT) NOPT=NMLQ
121: !120: !
122: !  QTEST is set for the final quenches with tighter convergence criteria.121: !  QTEST is set for the final quenches with tighter convergence criteria.
123: !122: !
124:       IF (QTEST) THEN123:       IF (QTEST) THEN
125:          GMAX=CQMAX124:          GMAX=CQMAX
126:       ELSE125:       ELSE
127:          GMAX=BQMAX126:          GMAX=BQMAX
128:       ENDIF127:       ENDIF
129: 128: 


r32722/takestep.f 2017-06-05 15:30:21.735812558 +0100 r32721/takestep.f 2017-06-05 15:30:23.855839914 +0100
 53: !        CLOSE(177) 53: !        CLOSE(177)
 54:  54: 
 55:       IF ((PERMOPT.OR.PERMINVOPT.OR.DISTOPT).AND.PERIODIC) RETURN  55:       IF ((PERMOPT.OR.PERMINVOPT.OR.DISTOPT).AND.PERIODIC) RETURN 
 56:       IF ((PERMOPT.OR.PERMINVOPT).AND.MKTRAPT) RETURN  56:       IF ((PERMOPT.OR.PERMINVOPT).AND.MKTRAPT) RETURN 
 57:        57:       
 58:       IF (MODEL1T) THEN 58:       IF (MODEL1T) THEN
 59:          RANDOM=(DPRAND()-0.5D0)*2.0D0 59:          RANDOM=(DPRAND()-0.5D0)*2.0D0
 60:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM 60:          COORDS(1,NP)=COORDSO(1,NP)+STEP(NP)*RANDOM
 61:          RETURN 61:          RETURN
 62:       ENDIF 62:       ENDIF
 63:       IF (MLQT.OR.MLPVB3T.OR.MLP3T.OR.ORBITALS) THEN 63:       IF (MLQT.OR.MLPVB3T.OR.MLP3T) THEN
 64:          DO J1=1,NATOMS 64:          DO J1=1,NATOMS
 65:             IF (FROZEN(J1)) CYCLE 65:             IF (FROZEN(J1)) CYCLE
 66:             RANDOM=(DPRAND()-0.5D0)*2.0D0 66:             RANDOM=(DPRAND()-0.5D0)*2.0D0
 67:             COORDS(J1,NP)=COORDS(J1,NP)+STEP(NP)*RANDOM 67:             COORDS(J1,NP)=COORDS(J1,NP)+STEP(NP)*RANDOM
 68:          ENDDO 68:          ENDDO
  69: !        WRITE(MYUNIT,'(A,F20.10)') 'coords(5)=',COORDS(5,1)
 69:          RETURN 70:          RETURN
 70:       ENDIF 71:       ENDIF
 71: ! 72: !
 72: !  Calling CENTRE if NORESET is .TRUE. can lead to problems with COORDSO containing an atom 73: !  Calling CENTRE if NORESET is .TRUE. can lead to problems with COORDSO containing an atom
 73: !  outside the permitted radius. Then it may be impossible to take a step that keeps all the 74: !  outside the permitted radius. Then it may be impossible to take a step that keeps all the
 74: !  atoms inside. 75: !  atoms inside.
 75: ! 76: !
 76:       PISQ = PI*PI 77:       PISQ = PI*PI
 77:       NTRIESMAX=100 78:       NTRIESMAX=100
 78:  79: 
 79: !     IF (CENT.AND.(.NOT.SEEDT)) CALL CENTRE2(COORDS(1:3*NATOMS,NP)) ! COORDS might have been shifted by symmetry 80: !     IF (CENT.AND.(.NOT.SEEDT)) CALL CENTRE2(COORDS(1:3*NATOMS,NP)) ! COORDS might have been shifted by symmetry
 80:       IF ((.NOT.NORESET).AND.(.NOT.PERMOPT).AND.(.NOT.DIFFRACTT).AND.(.NOT.BLNT).AND.(.NOT.PERIODIC)  81:       IF ((.NOT.NORESET).AND.(.NOT.PERMOPT).AND.(.NOT.DIFFRACTT).AND.(.NOT.BLNT).AND.(.NOT.PERIODIC) 
 81:      &   .AND.(.NOT.PERMINVOPT).AND.(.NOT.QCIPOTT).AND.(.NOT.MLP3T).AND.(.NOT.MLQT).AND.(.NOT.MKTRAPT).AND.(.NOT.MLPB3T) 82:      &   .AND.(.NOT.PERMINVOPT).AND.(.NOT.QCIPOTT).AND.(.NOT.MLP3T).AND.(.NOT.MLQT).AND.(.NOT.MKTRAPT).AND.(.NOT.MLPB3T)
 82:      &   .AND.(.NOT.GAUSST).AND.(.NOT.(CSMT.AND.(.NOT.SYMMETRIZECSM))).AND.(.NOT.PERCOLATET).AND.(.NOT.ORBITALS) 83:      &   .AND.(.NOT.GAUSST).AND.(.NOT.(CSMT.AND.(.NOT.SYMMETRIZECSM))).AND.(.NOT.PERCOLATET)
 83:      &   .AND.(.NOT.MLPVB3T)) THEN 84:      &   .AND.(.NOT.MLPVB3T)) THEN
 84: ! 85: !
 85: !        csw34> CHECK NOTHING HAS MOVED OUTSIDE THE CONTAINER RADIUS  86: !        csw34> CHECK NOTHING HAS MOVED OUTSIDE THE CONTAINER RADIUS 
 86: ! 87: !
 87:          DO J1=1,NATOMS 88:          DO J1=1,NATOMS
 88:             IF ((.NOT.(RIGID.OR.RIGIDINIT)).OR.(J1.LE.NATOMS/2)) THEN 89:             IF ((.NOT.(RIGID.OR.RIGIDINIT)).OR.(J1.LE.NATOMS/2)) THEN
 89:                J2=3*J1 90:                J2=3*J1
 90:                DUMMY2=COORDS(J2-2,NP)**2+COORDS(J2-1,NP)**2+COORDS(J2,NP)**2 91:                DUMMY2=COORDS(J2-2,NP)**2+COORDS(J2-1,NP)**2+COORDS(J2,NP)**2
 91:                IF (DUMMY2.GT.RADIUS) THEN 92:                IF (DUMMY2.GT.RADIUS) THEN
 92:                   IF (AMBERT) THEN ! jmc49 We don't really want a container at all in amber9, but this bit of code is being used  93:                   IF (AMBERT) THEN ! jmc49 We don't really want a container at all in amber9, but this bit of code is being used 
405:             IF (TMOVE(NP)) LOCALSTEP=STEP(NP)406:             IF (TMOVE(NP)) LOCALSTEP=STEP(NP)
406:             IF (TMOVE(NP).AND.TBP) LOCALSTEP=STEP(NP)*TBPINCR ! jdf43>407:             IF (TMOVE(NP).AND.TBP) LOCALSTEP=STEP(NP)*TBPINCR ! jdf43>
407:          ENDIF408:          ENDIF
408: !409: !
409: !  Angular move block.410: !  Angular move block.
410: !  If NORESET is .TRUE. then VAT won;t be set, so we should skip this block.411: !  If NORESET is .TRUE. then VAT won;t be set, so we should skip this block.
411: !412: !
412: !        IF (J1.EQ.JMAX) WRITE(MYUNIT,'(A,I6,4F15.5)') 'JMAX,VAT,ASTEP(NP),VMIN,prod=',JMAX,VAT(J1,NP), 413: !        IF (J1.EQ.JMAX) WRITE(MYUNIT,'(A,I6,4F15.5)') 'JMAX,VAT,ASTEP(NP),VMIN,prod=',JMAX,VAT(J1,NP), 
413: !    &                                    ASTEP(NP),VMIN,ASTEP(NP)*VMIN414: !    &                                    ASTEP(NP),VMIN,ASTEP(NP)*VMIN
414:          IF (((VAT(J1,NP).GT.ASTEP(NP)*VMIN).AND.(J1.EQ.JMAX)).AND.(.NOT.BLNT).AND.!(.NOT.RIGID).AND.415:          IF (((VAT(J1,NP).GT.ASTEP(NP)*VMIN).AND.(J1.EQ.JMAX)).AND.(.NOT.BLNT).AND.!(.NOT.RIGID).AND.
415:      &         (.NOT.DIFFRACTT).AND.(.NOT.GAUSST).AND.(.NOT.PERCOLATET).AND.(.NOT.MLPVB3T).AND.(.NOT.ORBITALS)416:      &         (.NOT.DIFFRACTT).AND.(.NOT.GAUSST).AND.(.NOT.PERCOLATET).AND.(.NOT.MLPVB3T)
416:      &        .AND.(.NOT.NORESET).AND.(.NOT.PERIODIC).AND.(.NOT.THOMSONT).AND.(.NOT.ONEDAPBCT).AND.(.NOT.ONEDPBCT)417:      &        .AND.(.NOT.NORESET).AND.(.NOT.PERIODIC).AND.(.NOT.THOMSONT).AND.(.NOT.ONEDAPBCT).AND.(.NOT.ONEDPBCT)
417:      &      .AND.(.NOT.TWODPBCT).AND.(.NOT.THREEDAPBCT).AND.(.NOT.THREEDPBCT).AND.(.NOT.QCIPOTT).AND.(.NOT.MLP3T).AND.(.NOT.MLPB3T)418:      &      .AND.(.NOT.TWODPBCT).AND.(.NOT.THREEDAPBCT).AND.(.NOT.THREEDPBCT).AND.(.NOT.QCIPOTT).AND.(.NOT.MLP3T).AND.(.NOT.MLPB3T)
418:      &       .AND.(.NOT.MLQT).AND.(.NOT.TWODAPBCT).AND.(.NOT.((NCORE(NP).GT.0).AND.(J1.GT.NATOMS-NCORE(NP))))) THEN419:      &       .AND.(.NOT.MLQT).AND.(.NOT.TWODAPBCT).AND.(.NOT.((NCORE(NP).GT.0).AND.(J1.GT.NATOMS-NCORE(NP))))) THEN
419: 420: 
420:             IF (DEBUG) WRITE(MYUNIT,'(A,I4,A,F12.4,A,F12.4,A,I4,A,F12.4)') 'angular move for atom ',J1, 421:             IF (DEBUG) WRITE(MYUNIT,'(A,I4,A,F12.4,A,F12.4,A,I4,A,F12.4)') 'angular move for atom ',J1, 
421:      &           ' V=',VMAX,' Vmin=',VMIN,' next most weakly bound atom is ',JMAX2,' V=',VMAX2422:      &           ' V=',VMAX,' Vmin=',VMIN,' next most weakly bound atom is ',JMAX2,' V=',VMAX2
422: 423: 
423:            THETA=DPRAND()*PI424:            THETA=DPRAND()*PI
424:            PHI=DPRAND()*PI*2.0D0425:            PHI=DPRAND()*PI*2.0D0
425: !426: !
599: !             COORDS(J2-1,NP)=COORDS(J2-1,NP)+LOCALSTEP*RANDOM*CMDIST(J1)/CMMAX600: !             COORDS(J2-1,NP)=COORDS(J2-1,NP)+LOCALSTEP*RANDOM*CMDIST(J1)/CMMAX
600:               COORDS(J2-1,NP)=COORDS(J2-1,NP)+LOCALSTEP*RANDOM*DUMMY601:               COORDS(J2-1,NP)=COORDS(J2-1,NP)+LOCALSTEP*RANDOM*DUMMY
601:               RANDOM=(DPRAND()-0.5D0)*2.0D0602:               RANDOM=(DPRAND()-0.5D0)*2.0D0
602: !             COORDS(J2,NP)=COORDS(J2,NP)+LOCALSTEP*RANDOM*CMDIST(J1)/CMMAX603: !             COORDS(J2,NP)=COORDS(J2,NP)+LOCALSTEP*RANDOM*CMDIST(J1)/CMMAX
603:               IF (.NOT.TWOD) COORDS(J2,NP)=COORDS(J2,NP)+LOCALSTEP*RANDOM*DUMMY604:               IF (.NOT.TWOD) COORDS(J2,NP)=COORDS(J2,NP)+LOCALSTEP*RANDOM*DUMMY
604:            ENDIF605:            ENDIF
605: !606: !
606: ! Stop atoms leaving the container in this step607: ! Stop atoms leaving the container in this step
607: !608: !
608:            IF ((.NOT.PERIODIC).AND.(.NOT.AMBERT).AND.(.NOT.(RIGID.AND.((J1.GT.NATOMS/2)))).AND.(.NOT.BLNT).AND.(.NOT.MLP3T)609:            IF ((.NOT.PERIODIC).AND.(.NOT.AMBERT).AND.(.NOT.(RIGID.AND.((J1.GT.NATOMS/2)))).AND.(.NOT.BLNT).AND.(.NOT.MLP3T)
609:      1     .AND.(.NOT.PERCOLATET).AND.(.NOT.DIFFRACTT).AND.(.NOT.THOMSONT).AND.(.NOT.GAUSST).AND.(.NOT.QCIPOTT).AND.(.NOT.ORBITALS)610:      1     .AND.(.NOT.PERCOLATET).AND.(.NOT.DIFFRACTT).AND.(.NOT.THOMSONT).AND.(.NOT.GAUSST).AND.(.NOT.QCIPOTT)
610:      2     .AND.(.NOT.MKTRAPT).AND.(.NOT.MLPB3T).AND.(.NOT.MLQT).AND.(.NOT.MLPVB3T)) THEN611:      2     .AND.(.NOT.MKTRAPT).AND.(.NOT.MLPB3T).AND.(.NOT.MLQT).AND.(.NOT.MLPVB3T)) THEN
611: !          IF ((.NOT.PERIODIC).AND.(.NOT.AMBER).AND.(.NOT.(RIGID.AND.(LOCALSTEP.EQ.0.0D0))).AND.(.NOT.BLNT)) THEN612: !          IF ((.NOT.PERIODIC).AND.(.NOT.AMBER).AND.(.NOT.(RIGID.AND.(LOCALSTEP.EQ.0.0D0))).AND.(.NOT.BLNT)) THEN
612:               DUMMY=COORDS(J2-2,NP)**2+COORDS(J2-1,NP)**2+COORDS(J2,NP)**2613:               DUMMY=COORDS(J2-2,NP)**2+COORDS(J2-1,NP)**2+COORDS(J2,NP)**2
613: !614: !
614: !  Simply rescaling the radius of an atom that leaves the container will bias the sampling615: !  Simply rescaling the radius of an atom that leaves the container will bias the sampling
615: !  of configuration space. However, we are not using takestep for bspt thermodynamic sampling!616: !  of configuration space. However, we are not using takestep for bspt thermodynamic sampling!
616: !  So, put the atom back in the container on the other side!617: !  So, put the atom back in the container on the other side!
617: !618: !
618: !              IF (DUMMY.GT.RADIUS) THEN619: !              IF (DUMMY.GT.RADIUS) THEN
619: !                 COORDS(J2-2,NP)=(SQRT(RADIUS)-0.5D0)*COORDS(J2-2,NP)/SQRT(DUMMY)620: !                 COORDS(J2-2,NP)=(SQRT(RADIUS)-0.5D0)*COORDS(J2-2,NP)/SQRT(DUMMY)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0