hdiff output

r31909/gay-berne.f90 2017-02-14 09:30:32.039027930 +0000 r31908/gay-berne.f90 2017-02-14 09:30:33.091041864 +0000
797: !        ELSE797: !        ELSE
798:          ENERGY = PYEPSNOT*ENERGY798:          ENERGY = PYEPSNOT*ENERGY
799:          G      = PYEPSNOT*G799:          G      = PYEPSNOT*G
800: 800: 
801: !        END IF801: !        END IF
802: !998     CONTINUE802: !998     CONTINUE
803:         RETURN803:         RETURN
804:       END SUBROUTINE PYGPERIODIC 804:       END SUBROUTINE PYGPERIODIC 
805: 805: 
806: 806: 
807:  
808: ! PY potential, dc430's implementation 
809: ! with PBC and continuous cutoff added 
810: ! plus extra LJ site 
811: ! 
812: ! This version is addressable with interactions scaled by the input epsilon matrix PYADDEPS 
813: ! 
814:       SUBROUTINE PYGPERIODICADD (X, G, ENERGY, GTEST, STEST) 
815:  
816:        use key, only: PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PCUTOFF,PARAMONOVCUTOFF,& 
817:                 &       pya1bin,pya2bin,pysignot,pyepsnot,radift,LJSITE,BLJSITE,PEPSILON1,& 
818:                 &       PSCALEFAC1,PSCALEFAC2,MAXINTERACTIONS,PYBINARYT,PYBINARYTYPE1,VT,PEPSILONATTR,PSIGMAATTR,PYADDEPS, NADDTARGET 
819:        use commons, only: natoms 
820:        use pymodule 
821:       IMPLICIT NONE 
822:  
823:       DOUBLE PRECISION :: X(3*NATOMS), G(3*NATOMS) 
824:       DOUBLE PRECISION :: RI(3), RJ(3), RIJ(3), NR(3), RIJSQ, ABSRIJ, P(3), THETA, THETA2, CT, ST 
825:       DOUBLE PRECISION :: AE1(3,3), BE1(3,3), AE2(3,3), BE2(3,3), RM(3,3), APB(3,3), APBINV(3,3) 
826:       DOUBLE PRECISION :: FCNT1, FCNT2, SRTFI1, SRTFI2, SRTFI(2), FMIN, LAMDAC1, LAMDAC2, ENERGY 
827:       DOUBLE PRECISION :: RHO1, RHO1SQ, RHO16, RHO112, RHO2, RHO2SQ, RHO26 
828:       DOUBLE PRECISION :: FCTR1, FCTR2, DVDF1, DVDF2  
829:       DOUBLE PRECISION :: DF1PI1, DF1PI2, DF1PI3, DF2PI1, DF2PI2, DF2PI3 
830:       DOUBLE PRECISION :: DF1PJ1, DF1PJ2, DF1PJ3, DF2PJ1, DF2PJ2, DF2PJ3  
831:       DOUBLE PRECISION :: RMI(3,3), RMJ(3,3), E(3,3), ESQ(3,3), DE1(3,3), DE2(3,3), DE3(3,3) 
832:       DOUBLE PRECISION :: DPI1RM(3,3), DPI2RM(3,3), DPI3RM(3,3), DPJ1RM(3,3), DPJ2RM(3,3), DPJ3RM(3,3) 
833:       DOUBLE PRECISION :: DF1DR(3), DF2DR(3), DG1DR(3), DG2DR(3) 
834:       DOUBLE PRECISION :: ARIBRJ(3), XC(3), XCMRI(3), XCMRJ(3), FIJ(3), TIJ(3), TJI(3) 
835:       double precision :: D1ABEZ(3,3), D2ABEZ(3,3), D3ABEZ(3,3), D1ABE(3,3), D2ABE(3,3), D3ABE(3,3)  
836:       LOGICAL          :: GTEST 
837:  
838:  
839: ! sf344 additions 
840:  
841:       DOUBLE PRECISION :: CLJ1(2),CLJ2(2),CLJ3(2),CLJ4(2),CLJ5(2),CLJ6(2),CLJ7(2),CLJ8(2),CLJ11(2),CLJ12(2), & 
842:                           & CLJ13(2),CLJ14(2),DFDR(2,3),DFP(2,6),dr(6),dCLJ1(2,12),dVDUMMY(12),LJ1(2),DUMMY,DUMMY1,DUMMY2,& 
843:                           & dLJ1(2,12),VDUMMY 
844:       DOUBLE PRECISION :: term2(maxinteractions),term3(maxinteractions), & 
845:                           & xlj(maxinteractions,2,3),rljvec(maxinteractions,3),rljunitvec(maxinteractions,3),& 
846:                           & drlj(maxinteractions,12),rlj(maxinteractions),rlj2(maxinteractions) 
847:       DOUBLE PRECISION :: LLJ(12,maxinteractions), dLLJ1(maxinteractions,12) 
848:       INTEGER          :: k, MJ1, MJ2 
849:       LOGICAL          :: STEST 
850:       DOUBLE PRECISION :: cosangle, sinangle, IMAT(3,3), tildematrix(3,3), rotmat(3,3), crossvector(3),tempcrd(6) 
851:      VT(1:NATOMS/2)=0.0D0 
852:      term2(:)=1.0D0 
853:      term3(:)=0.0D0 
854:  
855:          ENERGY = 0.D0 
856:          G(:)   = 0.D0 
857:          
858:         DO J1=1, REALNATOMS 
859:  
860:             J3      = 3*J1 
861:             J5      = OFFSET + J3 
862:             RI      = X(J3-2:J3) 
863:             P       = X(J5-2:J5) 
864:             angle=sqrt(dot_product(P,P)) 
865:             if(angle>twopi) then 
866: ! normalise angle-axis coordinates 
867:                 X(J5-2:J5)=X(J5-2:J5)/angle 
868:                 do 
869:                   angle=angle-twopi 
870:                   if(angle<2*pi) exit 
871:                 end do 
872: ! multiply with new angle 
873:                 X(J5-2:J5)=X(J5-2:J5)*angle 
874:             end if 
875:  
876:             CALL RMDRVT(P, RMIvec(J1,:,:), DPI1RMvec(J1,:,:), DPI2RMvec(J1,:,:), DPI3RMvec(J1,:,:), .TRUE.) 
877:  
878:         END DO 
879:  
880:          DO J1 = 1, REALNATOMS - 1 
881:  
882:             MJ1=MOD(J1-1,NADDTARGET)+1 
883:             J3      = 3*J1 
884:             J5      = OFFSET + J3 
885:             RI      = X(J3-2:J3) 
886:             P       = X(J5-2:J5) 
887: !     ROTATION MATRIX 
888:  
889: !            CALL RMDRVT(P, RMI, DPI1RM, DPI2RM, DPI3RM, GTEST) 
890:             RMI(:,:)=RMIvec(J1,:,:) 
891:             DPI1RM(:,:)=DPI1RMvec(J1,:,:) 
892:             DPI2RM(:,:)=DPI2RMvec(J1,:,:) 
893:             DPI3RM(:,:)=DPI3RMvec(J1,:,:) 
894:  
895:             AE1 = MATMUL(RMI,(MATMUL(AEZR1(J1,:,:),(TRANSPOSE(RMI))))) 
896:  
897:             IF (RADIFT) THEN 
898:  
899:                AE2 = MATMUL(RMI,(MATMUL(AEZR2(J1,:,:),(TRANSPOSE(RMI)))))          
900:  
901:             ENDIF 
902:  
903: !     BEGIN INNER LOOP OVER PARTICLES 
904:  
905:             DO J2 = J1 + 1, REALNATOMS 
906:  
907:                MJ2=MOD(J2-1,NADDTARGET)+1 
908:                J4     = 3*J2 
909:                J6     = OFFSET + J4 
910:                RJ     = X(J4-2:J4)  
911:                P      = X(J6-2:J6) 
912:  
913: !     ROTATION MATRIX 
914:  
915: !               CALL RMDRVT(P, RMJ, DPJ1RM, DPJ2RM, DPJ3RM, GTEST)                
916:                RMJ(:,:)=RMIvec(J2,:,:) 
917:                DPJ1RM(:,:)=DPI1RMvec(J2,:,:) 
918:                DPJ2RM(:,:)=DPI2RMvec(J2,:,:) 
919:                DPJ3RM(:,:)=DPI3RMvec(J2,:,:) 
920:       
921:                BE1 = MATMUL(RMJ,(MATMUL(AEZR1(J2,:,:),(TRANSPOSE(RMJ))))) 
922:  
923:                IF (RADIFT) THEN 
924:     
925:                   BE2 = MATMUL(RMJ,(MATMUL(AEZR2(J2,:,:),(TRANSPOSE(RMJ))))) 
926:  
927:                ENDIF 
928:  
929: !     CALCULATE SEPARATION 
930:  
931:                RIJ    = RI - RJ 
932:                RIJSQ  = DOT_PRODUCT(RIJ,RIJ) 
933:                ABSRIJ = DSQRT(RIJSQ) 
934:                NR     = RIJ / ABSRIJ 
935:  
936:                 IF(PARAMONOVCUTOFF.AND.RIJSQ>cut) GOTO 124 
937:  
938: !     CALCULATE ECF 
939:  
940:                CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE1, BE1, RIJ, LAMDAC1, FMIN) 
941:  
942:                FCNT1   = - FMIN 
943:                SRTFI1  = 1.D0 / DSQRT(FCNT1) 
944:                APB     = LAMDAC1 * AE1 + (1.D0 - LAMDAC1) * BE1 
945:                 
946:                CALL MTRXIN (APB, APBINV) 
947:  
948:                ARIBRJ =  LAMDAC1 * MATMUL(AE1,RI) + (1.D0 - LAMDAC1) * MATMUL(BE1,RJ) 
949:                XC     =  MATMUL(APBINV, ARIBRJ) 
950:                XCMRI  = XC - RI 
951:                XCMRJ  = XC - RJ 
952:                DF1DR  = - 2.D0 * LAMDAC1 * MATMUL(AE1,XCMRI) 
953:  
954:                D1ABEZ = MATMUL(DPI1RM,AEZR1(J1,:,:)) 
955:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
956:  
957:                D2ABEZ = MATMUL(DPI2RM,AEZR1(J1,:,:)) 
958:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
959:  
960:                D3ABEZ = MATMUL(DPI3RM,AEZR1(J1,:,:)) 
961:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
962:  
963:                DF1PI1 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D1ABE,XCMRI)) 
964:                DF1PI2 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D2ABE,XCMRI)) 
965:                DF1PI3 = LAMDAC1*DOT_PRODUCT(XCMRI,MATMUL(D3ABE,XCMRI)) 
966:  
967:                D1ABEZ = MATMUL(DPJ1RM,AEZR1(J2,:,:)) 
968:                D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
969:  
970:                D2ABEZ = MATMUL(DPJ2RM,AEZR1(J2,:,:)) 
971:                D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
972:  
973:                D3ABEZ = MATMUL(DPJ3RM,AEZR1(J2,:,:)) 
974:                D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ)))  
975:                 
976:                DF1PJ1 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D1ABE,XCMRJ)) 
977:                DF1PJ2 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D2ABE,XCMRJ)) 
978:                DF1PJ3 = (1.D0-LAMDAC1)*DOT_PRODUCT(XCMRJ,MATMUL(D3ABE,XCMRJ)) 
979:  
980:                RHO1   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI1 + PYSIGNOT) 
981:                RHO1SQ = RHO1*RHO1 
982:                RHO16  = RHO1SQ*RHO1SQ*RHO1SQ 
983:                RHO112 = RHO16 * RHO16 
984:  
985:                FCTR1  = 0.5D0*ABSRIJ*SRTFI1/(FCNT1*PYSIGNOT) 
986:                DG1DR  = (1.D0-SRTFI1)*NR/PYSIGNOT + FCTR1*DF1DR 
987:                DVDF1  = -2.D0*RHO112*RHO1*FCTR1 
988:  
989:                IF (RADIFT) THEN 
990:  
991:                   CALL BRENTMIN (0.D0, 0.51D0, 1.D0, AE2, BE2, RIJ, LAMDAC2, FMIN) 
992:  
993:                   FCNT2   = - FMIN 
994:                   SRTFI2  = 1.D0 / DSQRT(FCNT2) 
995:                   APB     = LAMDAC2 * AE2 + (1.D0 - LAMDAC2) * BE2 
996:  
997:                   CALL MTRXIN (APB, APBINV) 
998:  
999:                   ARIBRJ =  LAMDAC2 * MATMUL(AE2,RI) + (1.D0 - LAMDAC2) * MATMUL(BE2,RJ) 
1000:                   XC     =  MATMUL(APBINV, ARIBRJ) 
1001:                   XCMRI  = XC - RI 
1002:                   XCMRJ  = XC - RJ 
1003:                   DF2DR  = - 2.D0 * LAMDAC2 * MATMUL(AE2,XCMRI) 
1004:  
1005:                   RHO2   = PYSIGNOT / (ABSRIJ - ABSRIJ*SRTFI2 + PYSIGNOT) 
1006:                   RHO2SQ = RHO2*RHO2 
1007:                   RHO26  = RHO2SQ*RHO2SQ*RHO2SQ 
1008:                 
1009:                   FCTR2  = 0.5D0*ABSRIJ*SRTFI2/(FCNT2*PYSIGNOT) 
1010:                   DG2DR  = (1.D0-SRTFI2)*NR/PYSIGNOT+FCTR2*DF2DR 
1011:                   DVDF2  = RHO26*RHO2*FCTR2 
1012:  
1013:                   D1ABEZ = MATMUL(DPI1RM,AEZR2(J1,:,:)) 
1014:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D1ABEZ))) 
1015:  
1016:                   D2ABEZ = MATMUL(DPI2RM,AEZR2(J1,:,:)) 
1017:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D2ABEZ))) 
1018:  
1019:                   D3ABEZ = MATMUL(DPI3RM,AEZR2(J1,:,:)) 
1020:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMI))) + MATMUL(RMI,(TRANSPOSE(D3ABEZ))) 
1021:  
1022:                   DF2PI1 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D1ABE),XCMRI) 
1023:                   DF2PI2 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D2ABE),XCMRI) 
1024:                   DF2PI3 = LAMDAC2*DOT_PRODUCT(MATMUL(XCMRI,D3ABE),XCMRI) 
1025:  
1026:                   D1ABEZ = MATMUL(DPJ1RM,AEZR2(J2,:,:)) 
1027:                   D1ABE  = MATMUL(D1ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D1ABEZ))) 
1028:  
1029:                   D2ABEZ = MATMUL(DPJ2RM,AEZR2(J2,:,:)) 
1030:                   D2ABE  = MATMUL(D2ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D2ABEZ))) 
1031:  
1032:                   D3ABEZ = MATMUL(DPJ3RM,AEZR2(J2,:,:)) 
1033:                   D3ABE  = MATMUL(D3ABEZ,(TRANSPOSE(RMJ))) + MATMUL(RMJ,(TRANSPOSE(D3ABEZ))) 
1034:  
1035:                   DF2PJ1 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D1ABE),XCMRJ) 
1036:                   DF2PJ2 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D2ABE),XCMRJ) 
1037:                   DF2PJ3 = (1.D0-LAMDAC2)*DOT_PRODUCT(MATMUL(XCMRJ,D3ABE),XCMRJ) 
1038:  
1039:                ELSE 
1040:  
1041:                   RHO2   = RHO1 
1042:                   RHO26  = RHO16 
1043:                   DG2DR  = DG1DR 
1044:                   DVDF2  = RHO26*RHO2*FCTR1 
1045:                   DF2PI1 = DF1PI1 
1046:                   DF2PI2 = DF1PI2 
1047:                   DF2PI3 = DF1PI3 
1048:                   DF2PJ1 = DF1PJ1 
1049:                   DF2PJ2 = DF1PJ2 
1050:                   DF2PJ3 = DF1PJ3 
1051:  
1052:                   SRTFI2 = SRTFI1 
1053:                   DF2DR  = DF1DR 
1054:                   RHO2SQ = RHO1SQ 
1055:  
1056:                ENDIF 
1057:  
1058:  
1059: ! Correction terms to the potential if we require a cutoff at rc.  
1060:                  SRTFI(1)=SRTFI1 
1061:                  SRTFI(2)=SRTFI2 
1062:                  DFP(1,1)=DF1PI1 
1063:                  DFP(1,2)=DF1PI2 
1064:                  DFP(1,3)=DF1PI3 
1065:                  DFP(1,4)=DF1PJ1 
1066:                  DFP(1,5)=DF1PJ2 
1067:                  DFP(1,6)=DF1PJ3 
1068:                  DFP(2,1)=DF2PI1 
1069:                  DFP(2,2)=DF2PI2 
1070:                  DFP(2,3)=DF2PI3 
1071:                  DFP(2,4)=DF2PJ1 
1072:                  DFP(2,5)=DF2PJ2 
1073:                  DFP(2,6)=DF2PJ3 
1074:                  DFDR(1,:)=DF1DR(:) 
1075:                  DFDR(2,:)=DF2DR(:) 
1076:                  LJ1(1)=RHO1 
1077:                  LJ1(2)=RHO2 
1078:  !     CALCULATE PY POTENTIAL ENERGY 
1079: 123 CONTINUE  
1080:         VDUMMY = 0.0D0 
1081:                   
1082:                IF(LJSITE) THEN   
1083:                 do k=1,maxinteractions ! k=1 -- interaction between repulsive primary 'apex' sites 
1084:                          ! k=2 and k=3 -- interaction between secondary and primary 'apex' sites  
1085:                          ! k=4 -- interaction between secondary 'apex' sites (normal LJ interaction)  
1086:                         ! trying to modify code to allow for binary systems.  
1087:                         ! apex site heights will be defined in absolute units, 
1088:                         ! hence PYA1bin(J1,1) etc. will be removed from below 
1089:  
1090:                    IF(k==1) THEN 
1091:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1092:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1093:                    ELSE IF(k==2) THEN 
1094:                         DUMMY1=PSCALEFAC1vec(J1)!*PYA1bin(J1,1) 
1095:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1096:                    ELSE IF(k==3) THEN 
1097:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1098:                         DUMMY2=PSCALEFAC1vec(J2)!*PYA1bin(J2,1) 
1099:                    ELSE  
1100:                         DUMMY1=-PSCALEFAC2vec(J1)!*PYA1bin(J1,1) 
1101:                         DUMMY2=-PSCALEFAC2vec(J2)!*PYA1bin(J2,1) 
1102:                    END IF 
1103:                         ! first particle 
1104:                         xlj(k,1,:)=RI+DUMMY1*MATMUL(RMI,vecsbf)    ! vecsbf: (1,0,0) in the body frame of ellipsoid 
1105:  
1106:                         ! second particle 
1107:                         xlj(k,2,:)=RJ+DUMMY2*MATMUL(RMJ,vecsbf) 
1108:  
1109:                         ! separation between the LJ sites 
1110:                         rlj2(k)=(xlj(k,2,1)-xlj(k,1,1))**2+(xlj(k,2,2)-xlj(k,1,2))**2+(xlj(k,2,3)-xlj(k,1,3))**2 
1111:                         rlj(k)=sqrt(rlj2(k)) 
1112:                         rljvec(k,1)=xlj(k,2,1)-xlj(k,1,1) 
1113:                         rljvec(k,2)=xlj(k,2,2)-xlj(k,1,2) 
1114:                         rljvec(k,3)=xlj(k,2,3)-xlj(k,1,3) 
1115:  
1116:                         DUMMY=1.0D0/rlj(k) 
1117:                         rljunitvec(k,:)=rljvec(k,:)*DUMMY !/rlj(k) 
1118:  
1119:                         drlj(k,1)=DUMMY*(xlj(k,2,1)-xlj(k,1,1))         !drlj/dx1 
1120:                         drlj(k,2)=DUMMY*(xlj(k,2,2)-xlj(k,1,2))         !drlj/dy1 
1121:                         drlj(k,3)=DUMMY*(xlj(k,2,3)-xlj(k,1,3))         !drlj/dz1 
1122:                         drlj(k,4)=-drlj(k,1)                               !drlj/dx2 
1123:                         drlj(k,5)=-drlj(k,2)                               !drlj/dy2 
1124:                         drlj(k,6)=-drlj(k,3)                               !drlj/dz2 
1125:                         drlj(k,7) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI1RM,vecsbf)) !drlj/dpx1 
1126:                         drlj(k,8) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI2RM,vecsbf)) !drlj/dpy1 
1127:                         drlj(k,9) =-DUMMY*DUMMY1*DOT_PRODUCT(rljvec(k,:),MATMUL(DPI3RM,vecsbf)) !drlj/dpz1 
1128:                         drlj(k,10) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ1RM,vecsbf)) !drlj/dpx2 
1129:                         drlj(k,11) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ2RM,vecsbf)) !drlj/dpy2 
1130:                         drlj(k,12) =  DUMMY*DUMMY2*DOT_PRODUCT(rljvec(k,:),MATMUL(DPJ3RM,vecsbf)) !drlj/dpz2 
1131:  
1132:               ! interaction between the extra LJ sites: 
1133:                         LLJ(1,k)=sigma1(k)*DUMMY !/rlj(k) 
1134:                         LLJ(2,k)=LLJ(1,k)**2 
1135:                         LLJ(3,k)=LLJ(2,k)*LLJ(1,k) 
1136:                         LLJ(4,k)=LLJ(2,k)**2 
1137:                         LLJ(5,k)=LLJ(4,k)*LLJ(1,k) 
1138:                         LLJ(6,k)=LLJ(4,k)*LLJ(2,k) 
1139:                         LLJ(7,k)=LLJ(6,k)*LLJ(1,k) 
1140:                         LLJ(11,k)=LLJ(5,k)*LLJ(6,k) 
1141:                         LLJ(12,k)=LLJ(6,k)*LLJ(6,k) 
1142:  
1143:                             DO j=1,12 
1144:                                 dLLJ1(k,j) =-sigma1(k)*DUMMY*DUMMY*drlj(k,j) 
1145:                             END DO 
1146:  
1147:                  VDUMMY=VDUMMY+4.0D0*epsilon1(k,J1,J2)*term2(k)*(LLJ(12,k) - attr(k)*LLJ(6,k)) ! extra LJ sites 
1148:               end do ! k=1,4 
1149:              END IF ! IF(LJSITE) 
1150:                 IF (MJ1.EQ.MJ2) THEN ! same address - take just repulsion? 
1151:                    VDUMMY = VDUMMY + 4.0D0 * RHO112*PYADDEPS(MJ2,MJ1) 
1152:                ELSE 
1153:                    VDUMMY = VDUMMY + 4.0D0 * (RHO112 - 1.0D0 * RHO26)*PYADDEPS(MJ2,MJ1) 
1154:                ENDIF 
1155:  
1156:                ENERGY = ENERGY + VDUMMY 
1157:                VT(J1) = VT(J1) + VDUMMY 
1158:                VT(J2) = VT(J2) + VDUMMY        ! pair potentials 
1159:  
1160:             dVDUMMY(:) = 0.0D0 
1161:  
1162: !     CALCULATE GRADIENT 
1163:              FIJ = 0.0D0 
1164:              TIJ = 0.0D0 
1165:              TJI = 0.0D0 
1166:              dVDUMMY(:)=0.0D0 
1167:  
1168:            IF(LJSITE) THEN 
1169:             do k=1,maxinteractions 
1170:              do j=1,3 
1171:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1172:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently only repulsive 
1173:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1174:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1175:                FIJ(j) = dVDUMMY(j) 
1176:              end do 
1177:              do j=7,12 
1178:                dVDUMMY(j) = dVDUMMY(j) + 4.0D0*epsilon1(k,J1,J2)*(12.0D0*LLJ(11,k)*dLLJ1(k,j))*term2(k) + & 
1179:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(12,k)*term3(k)*drlj(k,j)! extra LJ site derivatives, currently only repulsive 
1180:                dVDUMMY(j) = dVDUMMY(j) - attr(k)*(4.0D0*epsilon1(k,J1,J2)*(6.0D0*LLJ(5,k)*dLLJ1(k,j))*term2(k) + & 
1181:                           & 4.0D0*epsilon1(k,J1,J2)*LLJ(6,k)*term3(k)*drlj(k,j)) ! attractive secondary apex site 
1182:              end do 
1183:             end do 
1184:            END IF 
1185:               IF (MJ1.EQ.MJ2) THEN ! same address - take just repulsion? 
1186:                    FIJ    = FIJ + 24.0D0*2.D0*RHO112*RHO1*DG1DR*PYADDEPS(MJ2,MJ1) 
1187:                    TIJ(1) = DVDF1*DF1PI1 
1188:                    TIJ(2) = DVDF1*DF1PI2 
1189:                    TIJ(3) = DVDF1*DF1PI3 
1190:                    TJI(1) = DVDF1*DF1PJ1 
1191:                    TJI(2) = DVDF1*DF1PJ2 
1192:                    TJI(3) = DVDF1*DF1PJ3 
1193:                    TIJ = dVDUMMY(7:9) + 24.0D0 * PYEPSNOT * TIJ*PYADDEPS(MJ2,MJ1) 
1194:                    TJI = dVDUMMY(10:12) + 24.0D0 * PYEPSNOT * TJI*PYADDEPS(MJ2,MJ1) 
1195:                ELSE 
1196:                    FIJ    = FIJ + 24.0D0*(2.D0*RHO112*RHO1*DG1DR - RHO26*RHO2*DG2DR)*PYADDEPS(MJ2,MJ1) 
1197:                    TIJ(1) = DVDF1*DF1PI1 + 1.0D0* DVDF2*DF2PI1 
1198:                    TIJ(2) = DVDF1*DF1PI2 + 1.0D0*  DVDF2*DF2PI2 
1199:                    TIJ(3) = DVDF1*DF1PI3 + 1.0D0*  DVDF2*DF2PI3 
1200:                    TJI(1) = DVDF1*DF1PJ1 + 1.0D0*  DVDF2*DF2PJ1 
1201:                    TJI(2) = DVDF1*DF1PJ2 + 1.0D0*  DVDF2*DF2PJ2 
1202:                    TJI(3) = DVDF1*DF1PJ3 + 1.0D0*  DVDF2*DF2PJ3 
1203:                    TIJ = dVDUMMY(7:9) + 24.0D0 * PYEPSNOT * TIJ*PYADDEPS(MJ2,MJ1) 
1204:                    TJI = dVDUMMY(10:12) + 24.0D0 * PYEPSNOT * TJI*PYADDEPS(MJ2,MJ1) 
1205:                ENDIF 
1206:  
1207:                G(J3-2:J3) = G(J3-2:J3) - FIJ 
1208:                G(J4-2:J4) = G(J4-2:J4) + FIJ 
1209:                G(J5-2:J5) = G(J5-2:J5) + TIJ 
1210:                G(J6-2:J6) = G(J6-2:J6) + TJI 
1211:  
1212: !     END INNER LOOP OVER PARTICLES 
1213: 124 CONTINUE 
1214:             ENDDO 
1215:  
1216: !     END OUTER LOOP OVER PARTICLES 
1217:  
1218:          ENDDO 
1219:          ENERGY = PYEPSNOT*ENERGY 
1220:          G      = PYEPSNOT*G 
1221:  
1222:         RETURN 
1223:       END SUBROUTINE PYGPERIODICADD 
1224:  
1225:  
1226: SUBROUTINE PARAMONOVNUMFIRSTDER(OLDX,STEST)807: SUBROUTINE PARAMONOVNUMFIRSTDER(OLDX,STEST)
1227: use commons, only : natoms808: use commons, only : natoms
1228: !use key809: !use key
1229: use modhess810: use modhess
1230: implicit none811: implicit none
1231: 812: 
1232: DOUBLE PRECISION  :: V(3*NATOMS),ENERGY,X(3*NATOMS),NUMGRAD(3*NATOMS),TRUEGRAD(3*NATOMS),&813: DOUBLE PRECISION  :: V(3*NATOMS),ENERGY,X(3*NATOMS),NUMGRAD(3*NATOMS),TRUEGRAD(3*NATOMS),&
1233:         &            OLDX(3*NATOMS),ETEMP(2,3*NATOMS),ksi814:         &            OLDX(3*NATOMS),ETEMP(2,3*NATOMS),ksi
1234: INTEGER    :: i,j815: INTEGER    :: i,j
1235: LOGICAL    :: GTEST,STEST816: LOGICAL    :: GTEST,STEST
1414:                         HESS(j,i)=HESS(i,j)995:                         HESS(j,i)=HESS(i,j)
1415:                 END DO996:                 END DO
1416: !        END IF997: !        END IF
1417: END DO998: END DO
1418: !WRITE(*,*) 'sf344> exiting Paramonovsecder', NATOMS, X(1)999: !WRITE(*,*) 'sf344> exiting Paramonovsecder', NATOMS, X(1)
1419: !WRITE(*,*) 'HESSIAN:'1000: !WRITE(*,*) 'HESSIAN:'
1420: !WRITE(*,'(12F10.3)') HESS(:,:)1001: !WRITE(*,'(12F10.3)') HESS(:,:)
1421: 1002: 
1422: END SUBROUTINE PYGPERIODICSECDER1003: END SUBROUTINE PYGPERIODICSECDER
1423: 1004: 
1424: SUBROUTINE PYGPERIODICSECDERADD(OLDX,STEST) 
1425: use commons 
1426: use modhess 
1427: implicit none 
1428:  
1429: DOUBLE PRECISION  :: V(3*NATOMS),EGB,X(3*NATOMS),OLDX(3*NATOMS),VTEMP(2,3*NATOMS),ksi 
1430: INTEGER    :: i,j 
1431: LOGICAL    :: GTEST,STEST 
1432:  
1433: ksi=0.00001D0 
1434: X(:)=OLDX(:) 
1435:  
1436: IF (.NOT.ALLOCATED(HESS)) ALLOCATE(HESS(3*NATOMS,3*NATOMS)) 
1437: HESS(:,:)=0.0D0 
1438: V(:)=0.0D0 
1439: VTEMP(:,:)=0.0D0 
1440:  
1441: DO i=1,3*NATOMS 
1442:  
1443:          X(i)=X(i)-ksi 
1444:   
1445:          CALL PYGPERIODICADD(X,V,EGB,GTEST,.FALSE.) 
1446:          VTEMP(1,:)=V(:) 
1447:   
1448:          X(i)=X(i)+2.0D0*ksi 
1449:  
1450:          CALL PYGPERIODICADD(X,V,EGB,GTEST,.FALSE.) 
1451:  
1452:          VTEMP(2,:)=V(:) 
1453:  
1454:                 DO j=i,3*NATOMS 
1455:                         HESS(i,j)=(VTEMP(2,j)-VTEMP(1,j))/(2.0D0*ksi) 
1456:                         HESS(j,i)=HESS(i,j) 
1457:                 END DO 
1458: END DO 
1459:  
1460: END SUBROUTINE PYGPERIODICSECDERADD 
1461:  
1462: SUBROUTINE EllipsoidsAAtoPolar(px1,py1,pz1,alpha,beta,gamma,alphadeg,betadeg,gammadeg)1005: SUBROUTINE EllipsoidsAAtoPolar(px1,py1,pz1,alpha,beta,gamma,alphadeg,betadeg,gammadeg)
1463: ! converts angle-axis coordinates px, py, pz to "polar-like" angles alpha and beta1006: ! converts angle-axis coordinates px, py, pz to "polar-like" angles alpha and beta
1464: ! px=cos(alpha)*cos(alpha)1007: ! px=cos(alpha)*cos(alpha)
1465: ! py=cos(alpha)*sin(beta)1008: ! py=cos(alpha)*sin(beta)
1466: ! pz=sin(alpha)1009: ! pz=sin(alpha)
1467: !use commons1010: !use commons
1468: implicit none1011: implicit none
1469: 1012: 
1470: DOUBLE PRECISION        :: px1, py1, pz1,px,py,pz1013: DOUBLE PRECISION        :: px1, py1, pz1,px,py,pz
1471: 1014: 


r31909/key.f90 2017-02-14 09:30:32.295031317 +0000 r31908/key.f90 2017-02-14 09:30:33.355045362 +0000
 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, & 49:      &        CONCUTFRACT, CONCUTABST, ENDNUMHESS2, CHARMMDFTBT, PAIRCOLOURT, REVERSEUPHILLT, WHOLEDNEB, &
 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, & 50:      &        NONEBMAX, READMASST, ONEDAPBCT, ONEDPBCT, INVTONEDPBCT, INVTTWODPBCT, TWODAPBCT, TWODPBCT, THREEDAPBCT, &
 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, & 51:      &        THREEDPBCT, FOURDAPBCT, FOURDPBCT, MODEDOWNT, CHEMSHIFT, TTM3T, &
 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, & 52:      &        NOINVERSION, INVERTPT, KNOWVECS, PMPATHT, AAORIENTT, MULTIJOBT, QUIPARGSTRT, QUIPPARAMST, HESSDUMPT, &
 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, & 53:      &        CLASSICALRATEST, TSPLITTINGT, HESSREADT, INSTANTONOPTT,INSTANTONSTARTDUMPT,VARSTEPOPTT, MOLPRO, REAXFFT, &
 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, & 54:      &        EIGENONLY,OVERCONV, GLJT,CLSTRINGT,CLSTRINGTST, PHI4MODT, EX1DT, MCPATHT, MCBIAST, RPHT, TWISTT, MCPATH2T, &
 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, & 55:      &        PBST, SSHT, GAUSSIAN03, CPPNEBT, CUDAT, CUDATIMET, TRUSTMODET,MODELOST, METRICTENSOR, INTSPRINGACTIVET, &
 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, & 56:      &        PERMGUESS, QCIPERMCHECK, DUMPFRQST, MULTIPOTT, MLP3T, MLPB3T, DUMPBESTPATH, ALIGNRBST, AVOID_COLLISIONS, MLPPROB, &
 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, & 57:      &        MALONALDEHYDE, MLPNEWREG, DJWRBT, STEALTHYT, STEALTV, LJADDT, MLPB3NEWT, MLPVB3T, &
 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, & 58:      &        QCIPOTT, QCIPOT2T, QCIRADSHIFTT, QCINOREPINT, QCIAMBERT, SLERPT, NOTRANSROTT, MAXGAPT, BULKBOXT, GDSQT, FLATTESTT, &
 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS, PYADDT 59:      &        MLQT, MLQPROB, LJADD2T, MACROIONT, NOREGBIAS
 60:  60: 
 61: ! sy349 > for testing the flatpath after dneb 61: ! sy349 > for testing the flatpath after dneb
 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:) 62:       !LOGICAL, ALLOCATABLE :: FLATPATHT(:)
 63:       LOGICAL FLATPATHT 63:       LOGICAL FLATPATHT
 64:  64: 
 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted) 65: ! bf269 > polymer in a pore (non-bonding (LJ) energy from neighbours is not subtracted)
 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential? 66:       LOGICAL :: PORE8T = .FALSE. ! add 8th power cylindrical pore to the potential?
 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z) 67:       INTEGER :: PORE8_AXIS = 3 ! principal axis of the cylindric pore (1:x, 2:y, 3:z)
 68:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1 68:       DOUBLE PRECISION :: PORE8_ENERGY = 1.0d1 ! energy of the pore when radius = 1
 69:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads 69:       LOGICAL :: HARMPOLYT = .FALSE. ! add harmonic bonds between the beads
107:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &107:      &        STOCKEXP, JPARAM, MCPATHTEMP, MCPATHDMAX, MCPATHSTEP, MCPATHACCRATIO, BIASFAC, &
108:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &108:      &        MCADDDEV, MCPATHQMIN, MCPATHQMAX, RPHQMIN, RPHQMAX, RPHTEMP, TWISTF, TWISTREF, MCPATHADDREF, &
109:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &109:      &        MCPATHGWS, MCPATHGWQ, MCPATHNEGLECT, MCPATHTOL, FRAMESDIFF,TMRATIO, INTMINFAC, MLPLAMBDA, COLL_TOL, KLIM, SCA, &
110:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2110:      &        NEBMAXERISE, GDSQ, FLATEDIFF, QCIADDREPCUT, QCIADDREPEPS, QCIRADSHIFT, INTCONCUT, MLQLAMBDA, FRQCONV, FRQCONV2
111: 111: 
112: !     sf344112: !     sf344
113:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &113:       DOUBLE PRECISION :: PCUTOFF,PYA11(3),PYA21(3),PYA12(3),PYA22(3),PEPSILON1(3),PSCALEFAC1(2),PSCALEFAC2(2), &
114:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &114:      &                     PEPSILONATTR(2),PSIGMAATTR(2), PYOVERLAPTHRESH, LJSITECOORDS(3), LJGSITESIGMA, LJGSITEEPS, &
115:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT115:      &                     PYLOCALSTEP(2),PYCFTHRESH,PYGRAVITYC1,PYGRAVITYC2,PERCCUT
116:  116:  
117:       DOUBLE PRECISION, ALLOCATABLE :: PYADDEPS(:,:) 
118:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)117:       DOUBLE PRECISION, ALLOCATABLE :: POINTSDECA(:), POINTSICOS(:)
119:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)118:       DOUBLE PRECISION, ALLOCATABLE :: VT(:), pya1bin(:,:),pya2bin(:,:),SITE(:,:)
120:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF119:       LOGICAL          :: LJSITE,BLJSITE,LJSITEATTR,PYBINARYT,PARAMONOVPBCX,PARAMONOVPBCY,PARAMONOVPBCZ,PARAMONOVCUTOFF
121:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION120:       LOGICAL          :: PYGPERIODICT,ELLIPSOIDT,LJSITECOORDST,REALIGNXYZ,MULTISITEPYT,LJGSITET,NORMALMODET,PYT,PYCOLDFUSION
122:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE121:       INTEGER          :: MAXINTERACTIONS,PYBINARYTYPE1,MYUNIT,NPYSITE
123:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR122:       DOUBLE PRECISION :: GBANISOTROPYR, GBWELLDEPTHR
124: !     CHRIS123: !     CHRIS
125: !     MCP124: !     MCP
126:       LOGICAL :: AMHT125:       LOGICAL :: AMHT
127:       CHARACTER(LEN=3) :: RES_TYPE_ALA='ALA'126:       CHARACTER(LEN=3) :: RES_TYPE_ALA='ALA'


r31909/keywords.f 2017-02-14 09:30:32.559034816 +0000 r31908/keywords.f 2017-02-14 09:30:33.619048857 +0000
729:          MIEF_PBCT=.FALSE.729:          MIEF_PBCT=.FALSE.
730:          MIEF_CUTT=.FALSE.730:          MIEF_CUTT=.FALSE.
731:          MIEF_BOX(1:3) = 1.0D9731:          MIEF_BOX(1:3) = 1.0D9
732:          MIEF_RCUT= 1.0D9732:          MIEF_RCUT= 1.0D9
733:          !733:          !
734:          ! UNDOCUMENTED keywords/parameters734:          ! UNDOCUMENTED keywords/parameters
735:          ! 735:          ! 
736:          TWISTT=.FALSE.736:          TWISTT=.FALSE.
737:          LJADDT=.FALSE.737:          LJADDT=.FALSE.
738:          LJADD2T=.FALSE.738:          LJADD2T=.FALSE.
739:          PYADDT=.FALSE. 
740:          INVERTPT=.FALSE.739:          INVERTPT=.FALSE.
741:          DNEBEFRAC=0.0D0740:          DNEBEFRAC=0.0D0
742:          MORPHT=.FALSE.741:          MORPHT=.FALSE.
743:          GREATCIRCLET=.FALSE.742:          GREATCIRCLET=.FALSE.
744:          MAXTSENERGY=1.0D100743:          MAXTSENERGY=1.0D100
745:          MAXBARRIER=1.0D100744:          MAXBARRIER=1.0D100
746:          MAXMAXBARRIER=1.0D100745:          MAXMAXBARRIER=1.0D100
747:          ReoptimiseEndpoints=.False.746:          ReoptimiseEndpoints=.False.
748:          ANGLEAXIS=.FALSE.747:          ANGLEAXIS=.FALSE.
749:          NFAILMAX=2748:          NFAILMAX=2
5560:            CALL READF(PYGRAVITYC2)5559:            CALL READF(PYGRAVITYC2)
5561:            EFIELDT=.TRUE.5560:            EFIELDT=.TRUE.
5562:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN5561:       ELSE IF (WORD.EQ.'PYOVERLAPTHRESH') THEN
5563:             CALL READF(PYOVERLAPTHRESH)5562:             CALL READF(PYOVERLAPTHRESH)
5564:             PYLOCALSTEP(:)=1.0D05563:             PYLOCALSTEP(:)=1.0D0
5565:             WRITE(*,'(A,F8.3)') 'keywords> ellipsoids considered to overlap for an ECF value of ', PYOVERLAPTHRESH5564:             WRITE(*,'(A,F8.3)') 'keywords> ellipsoids considered to overlap for an ECF value of ', PYOVERLAPTHRESH
5566:             IF(NITEMS.GT.2) THEN5565:             IF(NITEMS.GT.2) THEN
5567:                CALL READF(PYLOCALSTEP(1))5566:                CALL READF(PYLOCALSTEP(1))
5568:                CALL READF(PYLOCALSTEP(2))5567:                CALL READF(PYLOCALSTEP(2))
5569:             END IF5568:             END IF
5570:       ELSE IF (WORD.EQ.'PYADD') THEN 
5571:          PYADDT=.TRUE. 
5572:          CALL READI(NADDTARGET) 
5573:          WRITE(*,'(A,I6)') 'keyword> Target cluster size is ',NADDTARGET 
5574:          IF (MOD(NATOMS/2,NADDTARGET).NE.0) THEN 
5575:             WRITE(*,'(A,I6)') 'keyword> ERROR, target cluster size is not a factor of the number of PY particles ', 
5576:      &         NATOMS/2 
5577:             STOP 
5578:          ENDIF 
5579:          LUNIT=GETUNIT() 
5580:          OPEN(LUNIT,FILE='epsilon',STATUS='OLD') 
5581:          IF (.NOT.ALLOCATED(PYADDEPS)) ALLOCATE(PYADDEPS(NADDTARGET,NADDTARGET)) 
5582:          DO J1=1,NADDTARGET 
5583:             DO J2=1,NADDTARGET 
5584:                READ(LUNIT,*) PYADDEPS(J2,J1) 
5585:                WRITE(*,'(2I6,G20.10)') J1,J2,PYADDEPS(J2,J1) 
5586:             ENDDO 
5587:          ENDDO 
5588:          CLOSE(LUNIT) 
5589:       ELSE IF (WORD.EQ.'PYGPERIODIC') THEN5569:       ELSE IF (WORD.EQ.'PYGPERIODIC') THEN
5590:             NRBSITES = 25570:             NRBSITES = 2
5591:             ALLOCATE(RBSITE(NRBSITES,3))5571:             ALLOCATE(RBSITE(NRBSITES,3))
5592:             PYGPERIODICT = .TRUE.5572:             PYGPERIODICT = .TRUE.
5593: !            ANGLEAXIS2=.TRUE.5573: !            ANGLEAXIS2=.TRUE.
5594:             RBAAT=.TRUE.5574:             RBAAT=.TRUE.
5595:             CALL READF(PYA1(1))5575:             CALL READF(PYA1(1))
5596:             CALL READF(PYA1(2))5576:             CALL READF(PYA1(2))
5597:             CALL READF(PYA1(3))5577:             CALL READF(PYA1(3))
5598:             CALL READF(PYA2(1))5578:             CALL READF(PYA2(1))


r31909/potential.f 2017-02-14 09:30:32.835038472 +0000 r31908/potential.f 2017-02-14 09:30:33.887052409 +0000
1815:          ELSE IF (PYGT) THEN1815:          ELSE IF (PYGT) THEN
1816:             CALL PYG (COORDS,VNEW,ENERGY,GTEST,SSTEST)1816:             CALL PYG (COORDS,VNEW,ENERGY,GTEST,SSTEST)
1817:             IF (SSTEST) THEN1817:             IF (SSTEST) THEN
1818:                CALL PYGSECDER(COORDS,SSTEST)1818:                CALL PYGSECDER(COORDS,SSTEST)
1819:             END IF1819:             END IF
1820:             IF (PTEST) THEN1820:             IF (PTEST) THEN
1821:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1821:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1822:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1822:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1823:             END IF1823:             END IF
1824:          ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN1824:          ELSE IF (PYGPERIODICT.OR.PYBINARYT) THEN
1825:             IF (PYADDT) THEN1825:             CALL PYGPERIODIC (COORDS,VNEW,ENERGY,GTEST,SSTEST)
1826:                CALL PYGPERIODICADD(COORDS,VNEW,ENERGY,GTEST,SSTEST) 
1827:             ELSE 
1828:                CALL PYGPERIODIC (COORDS,VNEW,ENERGY,GTEST,SSTEST) 
1829:             ENDIF 
1830:             ! CALL PARAMONOVNUMFIRSTDER (COORDS,VNEW,ENERGY,GTEST,SSTEST)1826:             ! CALL PARAMONOVNUMFIRSTDER (COORDS,VNEW,ENERGY,GTEST,SSTEST)
1831:             ! STOP1827:             ! STOP
1832:             ! WRITE(*,*) 'VNEW='1828:             ! WRITE(*,*) 'VNEW='
1833:             ! DO J1=1,3*NATOMS1829:             ! DO J1=1,3*NATOMS
1834:             ! WRITE(*,*) VNEW(J1)1830:             ! WRITE(*,*) VNEW(J1)
1835:             ! END DO1831:             ! END DO
1836:             ! STOP1832:             ! STOP
1837:             IF (SSTEST) THEN1833:             IF (SSTEST) THEN
1838:                IF (PYADDT) THEN1834:                CALL PYGPERIODICSECDER(COORDS,SSTEST)
1839:                   CALL PYGPERIODICSECDERADD(COORDS,SSTEST) 
1840:                ELSE 
1841:                   CALL PYGPERIODICSECDER(COORDS,SSTEST) 
1842:                ENDIF 
1843:             END IF1835:             END IF
1844:             IF (PTEST) THEN1836:             IF (PTEST) THEN
1845:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '1837:                WRITE(*,10) ' potential> Energy for last cycle=',ENERGY,'         '
1846:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '1838:                WRITE(ESTRING,10) 'Energy for last cycle=',ENERGY,'         '
1847:             END IF1839:             END IF
1848:          ELSE IF (PYT) THEN1840:          ELSE IF (PYT) THEN
1849:           call py(COORDS,VNEW,ENERGY,.true.)1841:           call py(COORDS,VNEW,ENERGY,.true.)
1850:             IF (SSTEST) THEN1842:             IF (SSTEST) THEN
1851:                CALL py_secder(COORDS,SSTEST)1843:                CALL py_secder(COORDS,SSTEST)
1852:             END IF1844:             END IF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0