hdiff output

r32248/minpermdistrbcom.f90 2017-04-04 20:30:09.901856508 +0100 r32247/minpermdistrbcom.f90 2017-04-04 20:30:10.625865988 +0100
101:          CMY = CMY + DUMMYA(J2-1)101:          CMY = CMY + DUMMYA(J2-1)
102:          CMZ = CMZ + DUMMYA(J2)102:          CMZ = CMZ + DUMMYA(J2)
103:       ENDDO103:       ENDDO
104: 104: 
105:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS105:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS
106:       DO J1 = 1, NATOMS106:       DO J1 = 1, NATOMS
107:          J2 = 3*J1107:          J2 = 3*J1
108:          DUMMYA(J2-2) = DUMMYA(J2-2) - CMX108:          DUMMYA(J2-2) = DUMMYA(J2-2) - CMX
109:          DUMMYA(J2-1) = DUMMYA(J2-1) - CMY109:          DUMMYA(J2-1) = DUMMYA(J2-1) - CMY
110: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)110: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)
111:         IF(EFIELDT.AND.MULTISITEPYT) THEN 111: !         IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT) 
112:          DUMMYA(J2)   = DUMMYA(J2) 
113:         ELSE 
114:          DUMMYA(J2)   = DUMMYA(J2)   - CMZ112:          DUMMYA(J2)   = DUMMYA(J2)   - CMZ
115:         END IF 
116:       ENDDO113:       ENDDO
117: 114: 
118:       IF (EFIELDT .AND. INVERT ==-1) THEN115:       IF (EFIELDT .AND. INVERT ==-1) THEN
119:          RMATI(:,:) = 0.D0116:          RMATI(:,:) = 0.D0
120:          RMATI(1,1) = 1.D0; RMATI(2,2) =-1.D0; RMATI(3,3) = 1.D0117:          RMATI(1,1) = 1.D0; RMATI(2,2) =-1.D0; RMATI(3,3) = 1.D0
121:          DO J1 = 1, NATOMS118:          DO J1 = 1, NATOMS
122:             J2 = 3*J1119:             J2 = 3*J1
123:             DUMMYC(J2-2:J2) = MATMUL(RMATI,DUMMYA(J2-2:J2))120:             DUMMYC(J2-2:J2) = MATMUL(RMATI,DUMMYA(J2-2:J2))
124:          ENDDO121:          ENDDO
125:       ELSE122:       ELSE
444: 441: 
445:       IMPLICIT NONE442:       IMPLICIT NONE
446:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2443:       INTEGER          :: NATOMS, I, J, J1, J2, JMAX1, JMAX2, NORBIT1, NCHOOSE1, NORBIT2, NCHOOSE2
447:       DOUBLE PRECISION :: X(3*NATOMS), T1(3*NATOMS)444:       DOUBLE PRECISION :: X(3*NATOMS), T1(3*NATOMS)
448:       DOUBLE PRECISION :: XS(3*NATOMS), T1S(3*NATOMS), T2S(3*NATOMS), DIST(NATOMS)445:       DOUBLE PRECISION :: XS(3*NATOMS), T1S(3*NATOMS), T2S(3*NATOMS), DIST(NATOMS)
449:       DOUBLE PRECISION :: AX(3), P(3), Q2(4), ROTM(3,3), ROTMINV(3,3)446:       DOUBLE PRECISION :: AX(3), P(3), Q2(4), ROTM(3,3), ROTMINV(3,3)
450:       DOUBLE PRECISION :: THETA, THETAH, COST, SINT, COSTH, SINTH, ST, FCT447:       DOUBLE PRECISION :: THETA, THETAH, COST, SINT, COSTH, SINTH, ST, FCT
451:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, DUMMY, PROJ, DMAX2, CUTOFF1, DTEMP448:       DOUBLE PRECISION :: CMX, CMY, CMZ, DMAX, DUMMY, PROJ, DMAX2, CUTOFF1, DTEMP
452:       LOGICAL          :: DEBUG449:       LOGICAL          :: DEBUG
453: 450: 
 451: !      EFIELDT = .FALSE.
454:       CUTOFF1 = 1.D-03452:       CUTOFF1 = 1.D-03
455: !453: !
456: !     Move centre of mass to the origin.454: !     Move centre of mass to the origin.
457: !455: !
458:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0456:       CMX = 0.D0; CMY = 0.D0; CMZ = 0.D0
459:       DO I = 1, NATOMS457:       DO I = 1, NATOMS
460:          J = 3*I458:          J = 3*I
461:          CMX = CMX + X(J-2)459:          CMX = CMX + X(J-2)
462:          CMY = CMY + X(J-1)460:          CMY = CMY + X(J-1)
463:          CMZ = CMZ + X(J)461:          CMZ = CMZ + X(J)
464:       ENDDO462:       ENDDO
465: 463: 
466:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS464:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS
467:       DO I = 1, NATOMS465:       DO I = 1, NATOMS
468:          J = 3*I466:          J = 3*I
469:          XS(J-2) = X(J-2) - CMX467:          XS(J-2) = X(J-2) - CMX
470:          XS(J-1) = X(J-1) - CMY468:          XS(J-1) = X(J-1) - CMY
471: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)469: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)
472:         IF(EFIELDT.AND.MULTISITEPYT) THEN470:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          XS(J)   = X(J)   - CMZ
473:          XS(J)   = X(J)  
474:         ELSE 
475:          XS(J)   = X(J)   - CMZ 
476:         END IF 
477:       ENDDO471:       ENDDO
478: 472: 
479:       DMAX    = -1.D0473:       DMAX    = -1.D0
480:       NORBIT1 = 1474:       NORBIT1 = 1
481: 475: 
482:       IF (EFIELDT) THEN476:       IF (EFIELDT) THEN
483:          T1S(:) = XS(:)477:          T1S(:) = XS(:)
484:          GOTO 100478:          GOTO 100
485:       ENDIF479:       ENDIF
486: 480: 
602: 596: 
603:          ENDIF597:          ENDIF
604: 598: 
605:       ENDDO599:       ENDDO
606: 600: 
607: !601: !
608: !     and now rotate it into the xz plane.602: !     and now rotate it into the xz plane.
609: !603: !
610: 604: 
611:       CALL ROTATMXZ(NATOMS, JMAX2, T1S, DTEMP, DIST, Q2, .TRUE.)605:       CALL ROTATMXZ(NATOMS, JMAX2, T1S, DTEMP, DIST, Q2, .TRUE.)
 606: !      EFIELDT = .TRUE.
612: 607: 
613:       END SUBROUTINE ORIENTA608:       END SUBROUTINE ORIENTA
614: 609: 
615: !     ----------------------------------------------------------------------------------------------610: !     ----------------------------------------------------------------------------------------------
616: 611: 
617:       SUBROUTINE ROTATMXZ(NATOMS, JDO, T1S, PROJ, DIST, Q1, ROTT)612:       SUBROUTINE ROTATMXZ(NATOMS, JDO, T1S, PROJ, DIST, Q1, ROTT)
618: 613: 
619:       USE KEY, ONLY : EFIELDT614:       USE KEY, ONLY : EFIELDT
620: 615: 
621:       IMPLICIT NONE616:       IMPLICIT NONE
710:          CMX = CMX + T(J-2)705:          CMX = CMX + T(J-2)
711:          CMY = CMY + T(J-1)706:          CMY = CMY + T(J-1)
712:          CMZ = CMZ + T(J)707:          CMZ = CMZ + T(J)
713:       ENDDO708:       ENDDO
714:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS709:       CMX = CMX/NATOMS; CMY = CMY/NATOMS; CMZ = CMZ/NATOMS
715:       DO I = 1, NATOMS710:       DO I = 1, NATOMS
716:          J      = 3*I711:          J      = 3*I
717:          X(J-2) = T(J-2) - CMX712:          X(J-2) = T(J-2) - CMX
718:          X(J-1) = T(J-1) - CMY713:          X(J-1) = T(J-1) - CMY
719: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)714: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)
720:          IF(EFIELDT.AND.MULTISITEPYT)715:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          X(J)   = T(J)   - CMZ
721:           X(J)   = T(J) 
722:          ELSE 
723:           X(J)   = T(J)   - CMZ 
724:          END IF 
725:       ENDDO716:       ENDDO
726: 717: 
727: !     extract the rotation matrix corresponding to rotation via Q718: !     extract the rotation matrix corresponding to rotation via Q
728: 719: 
729:       CALL QROTMAT(Q2,RM)720:       CALL QROTMAT(Q2,RM)
730: 721: 
731:       DO I = 1, NATOMS722:       DO I = 1, NATOMS
732: 723: 
733:          J        = 3*I724:          J        = 3*I
734:          T1(1:3)  = MATMUL(RM(:,:), X(J-2:J))725:          T1(1:3)  = MATMUL(RM(:,:), X(J-2:J))
768:          CMXA = CMXA + XA(J2+1)759:          CMXA = CMXA + XA(J2+1)
769:          CMYA = CMYA + XA(J2+2)760:          CMYA = CMYA + XA(J2+2)
770:          CMZA = CMZA + XA(J2+3)761:          CMZA = CMZA + XA(J2+3)
771:       ENDDO762:       ENDDO
772:       CMXA = CMXA/NSIZE; CMYA = CMYA/NSIZE; CMZA = CMZA/NSIZE763:       CMXA = CMXA/NSIZE; CMYA = CMYA/NSIZE; CMZA = CMZA/NSIZE
773:       DO J1 = 1, NSIZE764:       DO J1 = 1, NSIZE
774:          J2 = 3*(J1-1)765:          J2 = 3*(J1-1)
775:          XA(J2+1) = XA(J2+1) - CMXA766:          XA(J2+1) = XA(J2+1) - CMXA
776:          XA(J2+2) = XA(J2+2) - CMYA767:          XA(J2+2) = XA(J2+2) - CMYA
777: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)768: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)
778:          IF(EFIELDT.AND.MULTISITEPYT)769:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          XA(J2+3) = XA(J2+3) - CMZA
779:           XA(J2+3) = XA(J2+3)  
780:          ELSE 
781:           XA(J2+3) = XA(J2+3) - CMZA 
782:          END IF 
783:       ENDDO770:       ENDDO
784: 771: 
785:       CMXB = 0.0D0; CMYB = 0.0D0; CMZB = 0.0D0772:       CMXB = 0.0D0; CMYB = 0.0D0; CMZB = 0.0D0
786:       DO J1 = 1, NSIZE773:       DO J1 = 1, NSIZE
787:          J2 = 3*(J1-1)774:          J2 = 3*(J1-1)
788:          CMXB = CMXB + XB(J2+1)775:          CMXB = CMXB + XB(J2+1)
789:          CMYB = CMYB + XB(J2+2)776:          CMYB = CMYB + XB(J2+2)
790:          CMZB = CMZB + XB(J2+3)777:          CMZB = CMZB + XB(J2+3)
791:       ENDDO778:       ENDDO
792:       CMXB = CMXB/NSIZE; CMYB = CMYB/NSIZE; CMZB = CMZB/NSIZE779:       CMXB = CMXB/NSIZE; CMYB = CMYB/NSIZE; CMZB = CMZB/NSIZE
793:       DO J1 = 1, NSIZE780:       DO J1 = 1, NSIZE
794:          J2 = 3*(J1-1)781:          J2 = 3*(J1-1)
795:          XB(J2+1) = XB(J2+1) - CMXB782:          XB(J2+1) = XB(J2+1) - CMXB
796:          XB(J2+2) = XB(J2+2) - CMYB783:          XB(J2+2) = XB(J2+2) - CMYB
797: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)784: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)
798:          IF(EFIELDT.AND.MULTISITEPYT)785:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          XB(J2+3) = XB(J2+3) - CMZB
799:           XB(J2+3) = XB(J2+3) 
800:          ELSE 
801:           XB(J2+3) = XB(J2+3) - CMZB 
802:          END IF 
803:       ENDDO786:       ENDDO
804: 787: 
805:       QMAT(1:4,1:4) = 0.0D0788:       QMAT(1:4,1:4) = 0.0D0
806: 789: 
807:       DO J1 = 1, NSIZE790:       DO J1 = 1, NSIZE
808:          J2 = 3*(J1-1)791:          J2 = 3*(J1-1)
809:          XM = XA(J2+1) - XB(J2+1)792:          XM = XA(J2+1) - XB(J2+1)
810:          YM = XA(J2+2) - XB(J2+2)793:          YM = XA(J2+2) - XB(J2+2)
811:          ZM = XA(J2+3) - XB(J2+3)794:          ZM = XA(J2+3) - XB(J2+3)
812:          XP = XA(J2+1) + XB(J2+1)795:          XP = XA(J2+1) + XB(J2+1)
905: 888: 
906:       ENDIF889:       ENDIF
907: 890: 
908:       DIST = SQRT(MINV)891:       DIST = SQRT(MINV)
909: 892: 
910:       DO J1 = 1, NATOMS893:       DO J1 = 1, NATOMS
911:          J2 = 3*(J1-1)894:          J2 = 3*(J1-1)
912:          RB(J2+1) = RB(J2+1) - CMXB895:          RB(J2+1) = RB(J2+1) - CMXB
913:          RB(J2+2) = RB(J2+2) - CMYB896:          RB(J2+2) = RB(J2+2) - CMYB
914: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)897: ! do not move along the z axis in case of a gravity field (currently implemented for MULTISITEPY)
915:          IF(EFIELDT.AND.MULTISITEPYT)898:          IF(.NOT.EFIELDT.AND..NOT.MULTISITEPYT)          RB(J2+3) = RB(J2+3) - CMZB
916:           RB(J2+3) = RB(J2+3) 
917:          ELSE 
918:           RB(J2+3) = RB(J2+3) - CMZB 
919:          END IF 
920:       ENDDO899:       ENDDO
921: 900: 
922:       CALL NEWROTGEOMA(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)901:       CALL NEWROTGEOMA(NATOMS,RB,Q2,RM,CMXA,CMYA,CMZA)
923: 902: 
924:       DEALLOCATE(XA,XB)903:       DEALLOCATE(XA,XB)
925: 904: 
926:       END SUBROUTINE MINDISTA905:       END SUBROUTINE MINDISTA
927: 906: 
928: !     ----------------------------------------------------------------------------------------------907: !     ----------------------------------------------------------------------------------------------
929: 908: 


r32248/rbperm.f90 2017-04-04 20:30:10.149859757 +0100 r32247/rbperm.f90 2017-04-04 20:30:10.853868977 +0100
 68:       COORDSCOMA(:) = COORDSA(1:3*NATOMS/2) 68:       COORDSCOMA(:) = COORDSA(1:3*NATOMS/2)
 69:       COORDSCOMB(:) = COORDSB(1:3*NATOMS/2) 69:       COORDSCOMB(:) = COORDSB(1:3*NATOMS/2)
 70:  70: 
 71:       IF (DEBUG) THEN 71:       IF (DEBUG) THEN
 72:          WRITE(*,'(3F20.10)')COORDSA(1:3*NATOMS) 72:          WRITE(*,'(3F20.10)')COORDSA(1:3*NATOMS)
 73:          WRITE(*,*) 73:          WRITE(*,*)
 74:          WRITE(*,'(3F20.10)')COORDSB(1:3*NATOMS) 74:          WRITE(*,'(3F20.10)')COORDSB(1:3*NATOMS)
 75:       ENDIF 75:       ENDIF
 76:  76: 
 77: !     do not call minpermdistrbCOM for multisite PY  77: !     do not call minpermdistrbCOM for multisite PY 
 78:       IF(MULTISITEPYT.OR.PYGPERIODICT) GOTO 100 78:       IF(MULTISITEPYT.OR.PYGPERIODICT.OR.SANDBOXT) GOTO 100
 79:         NATOMS = NATOMS/2 79:         NATOMS = NATOMS/2
 80:  80: 
 81:         CALL MINPERMDISTRBCOM(COORDSCOMB,COORDSCOMA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT) 81:         CALL MINPERMDISTRBCOM(COORDSCOMB,COORDSCOMA,DISTANCE,DIST2,QBEST,RMATBEST,DEBUG,BOXLX,BOXLY,BOXLZ,BULKT)
 82:  82: 
 83:         NATOMS = 2*NATOMS 83:         NATOMS = 2*NATOMS
 84:  84: 
 85:         IF (SQRT(DISTANCE) <= GEOMDIFFTOL) THEN 85:         IF (SQRT(DISTANCE) <= GEOMDIFFTOL) THEN
 86:          DISTANCE = SQRT(DISTANCE) 86:          DISTANCE = SQRT(DISTANCE)
 87:          IF (DEBUG) PRINT '(A)',' rbpermdist> minpermdistrbcom suggests identical' 87:          IF (DEBUG) PRINT '(A)',' rbpermdist> minpermdistrbcom suggests identical'
 88:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0 88:          CMAX = 0.0D0; CMAY = 0.0D0; CMAZ = 0.0D0
1406:       T1(1:3*NATOMS) = X(1:3*NATOMS)1406:       T1(1:3*NATOMS) = X(1:3*NATOMS)
1407: 1407: 
1408:       END SUBROUTINE RBROTXZ1408:       END SUBROUTINE RBROTXZ
1409: 1409: 
1410: !     ----------------------------------------------------------------------------------------------1410: !     ----------------------------------------------------------------------------------------------
1411: 1411: 
1412:       SUBROUTINE RBFRAME(C,XS,Q,RMATBEST)1412:       SUBROUTINE RBFRAME(C,XS,Q,RMATBEST)
1413: 1413: 
1414:       USE COMMONS, ONLY: NATOMS, NRBSITES1414:       USE COMMONS, ONLY: NATOMS, NRBSITES
1415:       USE KEY, ONLY    : DBPT, DBPTDT, DMBLPYT, NTIPT, MSSTOCKT, PAHAT, PATCHYDT, PAPT, PTSTSTT, NTSITES, NCARBON, &1415:       USE KEY, ONLY    : DBPT, DBPTDT, DMBLPYT, NTIPT, MSSTOCKT, PAHAT, PATCHYDT, PAPT, PTSTSTT, NTSITES, NCARBON, &
1416: &                        NPYSITE, MULTISITEPYT, PYT, SILANET, PYGPERIODICT, SANDBOXT !jwrm2> unused, STOCKAAT1416: &                        NPYSITE, MULTISITEPYT, PYT, SILANET, PYGPERIODICT !jwrm2> unused, STOCKAAT
1417:       USE pymodule, ONLY : SITECOORDS,ELLST1,ELLMAT1417:       USE pymodule, ONLY : SITECOORDS,ELLST1,ELLMAT
1418: 1418: 
1419:       IMPLICIT NONE1419:       IMPLICIT NONE
1420:       1420:       
1421:       INTEGER :: J1, J2, J31421:       INTEGER :: J1, J2, J3
1422:       DOUBLE PRECISION :: XS(3*NTSITES)1422:       DOUBLE PRECISION :: XS(3*NTSITES)
1423:       CHARACTER(LEN=132),INTENT(IN) :: C1423:       CHARACTER(LEN=132),INTENT(IN) :: C
1424:       DOUBLE PRECISION :: Q(3*NATOMS),RMATBEST(3,3)1424:       DOUBLE PRECISION :: Q(3*NATOMS),RMATBEST(3,3)
1425: 1425: 
1426:       IF (DBPTDT ) THEN1426:       IF (DBPTDT ) THEN
1565:                         'ellipse ',ELLST1(J3,1)*2.0D0,ELLST1(J3,2)*2.0D0,ELLST1(J3,3)*2.0D0,&1565:                         'ellipse ',ELLST1(J3,1)*2.0D0,ELLST1(J3,2)*2.0D0,ELLST1(J3,3)*2.0D0,&
1566:                         ELLMAT(J3,1,1),ELLMAT(J3,1,2),ELLMAT(J3,1,3),ELLMAT(J3,2,1),ELLMAT(J3,2,2),ELLMAT(J3,2,3),&1566:                         ELLMAT(J3,1,1),ELLMAT(J3,1,2),ELLMAT(J3,1,3),ELLMAT(J3,2,1),ELLMAT(J3,2,2),ELLMAT(J3,2,3),&
1567:                         ELLMAT(J3,3,1),ELLMAT(J3,3,2),ELLMAT(J3,3,3),&1567:                         ELLMAT(J3,3,1),ELLMAT(J3,3,2),ELLMAT(J3,3,3),&
1568:                         'atom_vector',Q(3*NATOMS/2+3*(J2-1)+1),&1568:                         'atom_vector',Q(3*NATOMS/2+3*(J2-1)+1),&
1569:                         Q(3*NATOMS/2+3*(J2-1)+2),Q(3*NATOMS/2+3*(J2-1)+3)1569:                         Q(3*NATOMS/2+3*(J2-1)+2),Q(3*NATOMS/2+3*(J2-1)+3)
1570:                     END DO1570:                     END DO
1571:              END DO1571:              END DO
1572:       ELSEIF (PYT) THEN ! sf344> implementation still buggy, use MULTISITEPY instead!1572:       ELSEIF (PYT) THEN ! sf344> implementation still buggy, use MULTISITEPY instead!
1573:          WRITE(55,'(i6)') NATOMS*NRBSITES/21573:          WRITE(55,'(i6)') NATOMS*NRBSITES/2
1574:          WRITE(55,'(1x,a)') trim(adjustl(c))1574:          WRITE(55,'(1x,a)') trim(adjustl(c))
1575: !        WRITE(*,*) 'in here'1575:         WRITE(*,*) 'in here'
1576:          CALL PY_OUTPUT(55,Q)1576:          CALL PY_OUTPUT(55,Q)
1577: !      ELSEIF (STOCKAAT) THEN1577: !      ELSEIF (STOCKAAT) THEN
1578: 1578: 
1579: !         WRITE(55,'(i6)') (NATOMS/2)*NRBSITES1579: !         WRITE(55,'(i6)') (NATOMS/2)*NRBSITES
1580: !         WRITE(55,'(1x,a)') trim(adjustl(c))1580: !         WRITE(55,'(1x,a)') trim(adjustl(c))
1581: !         DO J1 = 1, NATOMS/21581: !         DO J1 = 1, NATOMS/2
1582: !            J2 = (J1-1)*31582: !            J2 = (J1-1)*3
1583: !            WRITE(55,'(A5,1X,3F20.10)') 'O ', XS(J2+1), XS(J2+2), XS(J2+3)1583: !            WRITE(55,'(A5,1X,3F20.10)') 'O ', XS(J2+1), XS(J2+2), XS(J2+3)
1584: !         ENDDO1584: !         ENDDO
1585:       ELSEIF (SANDBOXT) THEN 
1586:          WRITE(55,'(i6)') NTSITES 
1587:          WRITE(55,'(1x,a)') trim(adjustl(c)) 
1588:          CALL SANDBOX_OUTPUT(55,XS) 
1589: 1585: 
1590:       ELSEIF (SILANET) THEN1586:       ELSEIF (SILANET) THEN
1591: 1587: 
1592:          WRITE(55,'(i6)') (NATOMS/2)*NRBSITES1588:          WRITE(55,'(i6)') (NATOMS/2)*NRBSITES
1593:          WRITE(55,'(1x,a)') trim(adjustl(c))1589:          WRITE(55,'(1x,a)') trim(adjustl(c))
1594:          DO J1 = 1, NATOMS/21590:          DO J1 = 1, NATOMS/2
1595:             DO J2 = 1, NRBSITES1591:             DO J2 = 1, NRBSITES
1596:                J3 = ((J1-1)*NRBSITES + (J2-1))*31592:                J3 = ((J1-1)*NRBSITES + (J2-1))*3
1597:                IF (J2 == 1) THEN1593:                IF (J2 == 1) THEN
1598:                   WRITE(55,'(A5,1X,3F20.10)') 'Si ', XS(J3+1), XS(J3+2), XS(J3+3)1594:                   WRITE(55,'(A5,1X,3F20.10)') 'Si ', XS(J3+1), XS(J3+2), XS(J3+3)


r32248/sandbox.f90 2017-04-04 20:30:10.381862796 +0100 r32247/sandbox.f90 2017-04-04 20:30:11.077871909 +0100
870: !        STOP870: !        STOP
871:         do moli_index = 1, num_mols871:         do moli_index = 1, num_mols
872:             moli => molecules(moli_index)872:             moli => molecules(moli_index)
873:             if(frozen(moli_index)) grad(moli%r_index: moli%r_index + 2) = 0.d0873:             if(frozen(moli_index)) grad(moli%r_index: moli%r_index + 2) = 0.d0
874:             if(frozen(moli_index + num_mols)) grad(moli%p_index: moli%p_index + 2) = 0.d0874:             if(frozen(moli_index + num_mols)) grad(moli%p_index: moli%p_index + 2) = 0.d0
875:         end do875:         end do
876:     end if876:     end if
877: 877: 
878: end subroutine sandbox878: end subroutine sandbox
879: 879: 
880: subroutine sandbox_output(sandout_unit,xs)880: subroutine sandbox_output
881:     use commons, only: natoms 881:     use commons, only: natoms 
882:     use sandbox_module882:     use sandbox_module
883:     implicit none883:     implicit none
884: 884: 
885:     integer :: sandout_unit, coords_unit,j1,j2885:     integer :: sandout_unit, coords_unit
886:     integer :: min_num, mol_num, site_num, atom_num, atom_index886:     integer :: min_num, mol_num, site_num, atom_num, atom_index
887:     type(sandbox_molecule), pointer :: mol887:     type(sandbox_molecule), pointer :: mol
888:     type(sandbox_site), pointer :: site888:     type(sandbox_site), pointer :: site
889:     character(len=25) :: sandout_name, coords_name, min_name, node889:     character(len=25) :: sandout_name, coords_name, min_name, node
890:     double precision :: rotmat(3, 3), dummy(3, 3),xs(3*num_atoms)890:     double precision :: rotmat(3, 3), dummy(3, 3)
891:     double precision :: qminp1(3*natoms) 
892: 891: 
893:     integer :: getunit892:     integer :: getunit
894: 893: 
895:       integer,allocatable :: FF(:),INTEFF(:) ! NSAVE894:       integer,allocatable :: FF(:),INTEFF(:) ! NSAVE
896:       INTEGER, ALLOCATABLE :: NPCALL_QMIN(:) ! NSAVE895:       INTEGER, ALLOCATABLE :: NPCALL_QMIN(:) ! NSAVE
897:       DOUBLE PRECISION, ALLOCATABLE :: QMIN(:), INTEQMIN(:) ! NSAVE896:       DOUBLE PRECISION, ALLOCATABLE :: QMIN(:), INTEQMIN(:) ! NSAVE
898:       DOUBLE PRECISION, ALLOCATABLE :: QMINP(:,:), INTEQMINP(:,:)897:       DOUBLE PRECISION, ALLOCATABLE :: QMINP(:,:), INTEQMINP(:,:)
899:       INTEGER, ALLOCATABLE :: QMINT(:,:), QMINNATOMS(:)898:       INTEGER, ALLOCATABLE :: QMINT(:,:), QMINNATOMS(:)
900: 899: 
901: 900: 
902: !    sandout_unit = getunit()901:     sandout_unit = getunit()
 902: 
 903:     ! open sandout file for writing
 904:     if(mpit) then
 905:         write(node, *) mynode + 1
 906:         sandout_name = 'sandout.' // trim(adjustl(node)) // '.xyz'
 907:     else
 908:         sandout_name = 'sandout.xyz'
 909:     end if
 910:     open(unit=sandout_unit, file=sandout_name, status='unknown')
 911: 
 912:     ! loop over saved minima
 913:     do min_num = 1, nsave
 914:         ! put number of atoms and comment line. for now this we'll do just one atom per molecule
 915:         !write(sandout_unit, '(I0)') num_mols
 916:         write(sandout_unit, '(I0)') num_atoms
 917: !        write(sandout_unit, '(A,1x,I0,A,1x,F20.10,1x,A,1x,I0)') 'energy of minimum', min_num, '=', qmin(min_num), 'first found at step', ff(min_num)
903: 918: 
904:         j1=0 
905:         ! loop over molecules919:         ! loop over molecules
906:         do mol_num = 1, size(molecules)920:         do mol_num = 1, size(molecules)
907:             mol => molecules(mol_num)921:             mol => molecules(mol_num)
908: 922: 
909: !            call rmdrvt(qminp1(mol%p_index: mol%p_index + 2), rotmat, dummy, dummy, dummy, .false.)923:             call rmdrvt(qminp(min_num, mol%p_index: mol%p_index + 2), rotmat, dummy, dummy, dummy, .false.)
910: 924: 
911: !            call update_sandbox_molecule(.false., mol, qminp1(mol%r_index: mol%r_index + 2),qminp1(mol%p_index: mol%p_index + 2))925:             call update_sandbox_molecule(.false., mol, qminp(min_num, mol%r_index: mol%r_index + 2),qminp(min_num,mol%p_index: mol%p_index + 2))
912: 926: 
913:             !write(sandout_unit, '(A6,3F10.5,A12,1X,3F10.5)') mol%atom_symbol, qminp(min_num, mol%r_index: mol%r_index + 2), &927:             !write(sandout_unit, '(A6,3F10.5,A12,1X,3F10.5)') mol%atom_symbol, qminp(min_num, mol%r_index: mol%r_index + 2), &
914:              !   & 'atom_vector', matmul(rotmat, z_hat)928:              !   & 'atom_vector', matmul(rotmat, z_hat)
915:             do atom_num=1,size(molecules(mol_num)%sites)929:             do atom_num=1,size(molecules(mol_num)%sites)
916:                j2=3*j1930:                write(sandout_unit, '(A6,3F10.5,A12,1X,3F10.5)')molecules(mol_num)%sites(atom_num)%name,molecules(mol_num)%sites(atom_num)%r
917:                write(sandout_unit, '(A6,3F10.5,A12,1X,3F10.5)')molecules(mol_num)%sites(atom_num)%name, XS(J2+1), XS(J2+2), XS(J2+3) 
918:                j1=j1+1 
919:             enddo931:             enddo
920:         end do932:         end do
921: !         DO J1 = 1, NATOMS/2 
922: !            J2 = (J1-1)*3 
923: !            WRITE(55,'(A5,1X,3F20.10)') 'O ', XS(J2+1), XS(J2+2), XS(J2+3) 
924: !         ENDDO 
925: 933: 
 934:         ! output coords. files
 935:         coords_unit = getunit()
 936:         write(min_name, '(i3)') min_num
 937: 
 938:         if(mpit) then
 939:             coords_name = 'coords.' // trim(adjustl(min_name)) // '.' // trim(adjustl(node))
 940:         else
 941:             coords_name = 'coords.' // trim(adjustl(min_name))
 942:         end if
 943: 
 944:         open(coords_unit, file=coords_name, status='unknown')
 945:         do atom_num = 1, natoms
 946:             atom_index = 3 * (atom_num - 1) + 1
 947:             write(coords_unit, *) qminp(min_num, atom_index), qminp(min_num, atom_index + 1), qminp(min_num, atom_index + 2)
 948:         end do
 949:         close(coords_unit)
 950: 
 951:     end do
926: 952: 
927: !    close(sandout_unit)953:     close(sandout_unit)
928:             954:             
929: end subroutine sandbox_output955: end subroutine sandbox_output
930: 956: 
931: subroutine sandbox_takestep(np)957: subroutine sandbox_takestep(np)
932:     use key, only: myunit,twod, frozen958:     use key, only: myunit,twod, frozen
933:     use sandbox_module959:     use sandbox_module
934:     use vec3, only: vec_len, vec_random960:     use vec3, only: vec_len, vec_random
935:     use rotations, only: rot_takestep_aa961:     use rotations, only: rot_takestep_aa
936: 962: 
937:     implicit none963:     implicit none


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0