PROGRAM CONST C C CREATE CONSTRAINTS MATRIX GIVEN SOME HINTS C C MUST COMPILE WITH REAL -> DOUBLE PROMOTION C C LAST REVISED 17/2/97 MJR C 26/3/97 SO THAT 1-1E-7 IS RECOGNISED AS 0 C IMPLICIT NONE INTEGER MXSPEC,MXIONS PARAMETER(MXSPEC=10,MXIONS=50) INTEGER NSPEC,NIONS(MXSPEC),RMOVE(MXSPEC,MXIONS),AMOVE(3) INTEGER I,J,ICHOICE,INVERT(MXSPEC,MXIONS),CAXIS(3) INTEGER NOFF(MXSPEC),LINE(3*MXSPEC*MXIONS) INTEGER NAT,DOF,NIONST,K,DOFPA,NSP,T,IATOM,NMOVE,NDOF LOGICAL CETEP,ZONLY REAL POSION(MXSPEC,MXIONS,3),CENTRE(3),RIMAGE(3) REAL HESS(3*MXIONS,3*MXIONS) REAL DUMMY,DIST2 C WRITE(*,*)'PLEASE INPUT NSPEC' READ(*,*)NSPEC IF (NSPEC.GT.MXSPEC) STOP 'NSPEC>MXSPEC' C WRITE(*,*)'PLEASE INPUT NIONS(I),I=1,NSPEC' READ(*,*)(NIONS(I),I=1,NSPEC) C WRITE(*,*)'PLEASE CHOOSE FROM:' WRITE(*,*)' 1) MOVE ALL ATOMS' WRITE(*,*)' 2) FIX JUST ONE ATOM' WRITE(*,*)' 3) CONSIDER RMOVES IN A FORT.15' C 100 READ(*,*)ICHOICE IF ((ICHOICE.LT.1).OR.(ICHOICE.GT.3)) GOTO 100 C DO I=1,MXSPEC DO J=1,MXIONS RMOVE(I,J)=0.0 INVERT(I,J)=-9999 END DO END DO C IF ((ICHOICE.EQ.1).OR.(ICHOICE.EQ.2)) THEN DO I=1,NSPEC DO J=1,NIONS(I) RMOVE(I,J)=1.0 END DO END DO IF (ICHOICE.EQ.2) THEN WRITE(*,*)'FIX WHICH ATOM (NSP NATOM)' READ(*,*) I,J RMOVE(I,J)=0 ENDIF ELSE OPEN(15,STATUS='OLD') REWIND(15) DO I=1,6 READ(15,*)DUMMY END DO DO I=1,NSPEC DO J=1,NIONS(I) READ(15,*)DUMMY,DUMMY,DUMMY,RMOVE(I,J) END DO DO J=1,NIONS(I) READ(15,*)DUMMY END DO END DO CLOSE(15) WRITE(*,*)'FORT.15 READ' ENDIF C WRITE(*,*)'IMPOSE INVERSION SYMMETRY? (1=YES, 0=NO)' READ(*,*)ICHOICE C IF(ICHOICE.EQ.1) THEN WRITE(*,*)'ABOUT WHICH POINT?' READ(*,*)(CENTRE(I),I=1,3) C OPEN(15,STATUS='OLD') REWIND(15) DO I=1,6 READ(15,*)DUMMY END DO DO I=1,NSPEC DO J=1,NIONS(I) READ(15,*)(POSION(I,J,K),K=1,3) END DO DO J=1,NIONS(I) READ(15,*)DUMMY END DO END DO CLOSE(15) C C GET ALL IONS IN 0<=X<1 RANGE C DO I=1,NSPEC DO J=1,NIONS(I) DO K=1,3 POSION(I,J,K)=MOD(POSION(I,J,K),1.0) IF (POSION(I,J,K).GT.1-1E-7) POSION(I,J,K)=0.0 ENDDO END DO END DO C C INVERT(NSPEC,NIONS) IS CONTRUCTED THUS: C C INVERT(NSP,I)=I IF NO INVERSION OR I MAPS TO I C INVERT(NSP,I)=J ) IF I MAPS TO J C INVERT(NSP,J)=-I ) C C THUS SCANNING INVERT(NSP,I).GT.0 GUARANTEES ALL PAIRS ARE HIT C ONCE ONLY C C DO I=1,NSPEC DO J=1,NIONS(I) IF(INVERT(I,J).EQ.-9999) THEN DO K=1,3 RIMAGE(K)= & MOD((-(POSION(I,J,K)-CENTRE(K))+CENTRE(K))+1.0,1.0) C FORCE WRAP AS REQUIRED IF (RIMAGE(K).GT.1-1E-7) RIMAGE(K)=0.0 END DO IATOM=0 DO T=1,NIONS(I) DIST2=0 DO K=1,3 DIST2=DIST2+(POSION(I,T,K)-RIMAGE(K))**2 END DO IF (DIST2.LE.1E-10) IATOM=T END DO IF (IATOM.EQ.0) THEN WRITE(*,*) & 'ERROR: CANNOT FIND IMAGE OF NSPEC=',I,' ATOM=',J WRITE(*,*)RIMAGE STOP ENDIF IF (RMOVE(I,J).NE.RMOVE(I,IATOM)) THEN WRITE(*,*)'WARNING: RMOVE NOT CONSITENT WITH INVERSION' WRITE(*,*)' FIXING ATOM' RMOVE(I,J)=0 RMOVE(I,IATOM)=0 ENDIF WRITE(*,*)i,j,' <--> ',i,iatom c write(*,1234)(posion(i,j,k),k=1,3), c & (posion(i,iatom,k),k=1,3) c1234 format(3F8.3,' <--> ',3F8.3) INVERT(I,IATOM)=-J INVERT(I,J)=IATOM ENDIF END DO END DO ELSE DO I=1,NSPEC DO J=1,NIONS(I) INVERT(I,J)=J END DO END DO END IF C WRITE(*,*)'PLEASE SPECIFY PERMITED MOTION DIRECTIONS' WRITE(*,*) & 'INPUT 1 OR 0 FOR EACH DIRECTION IN WHICH ATOMS SHOULD MOVE' WRITE(*,*)' EG 1 1 1 FOR FREE MOTION, 0 0 1 FOR MOTION IN Z ONLY' 111 READ(*,*) (AMOVE(I),I=1,3) C C NOW FOR SOME SANITY CHECKS C NMOVE=0 DOF=0 DOFPA=0 DO I=1,3 IF (AMOVE(I).NE.0) DOFPA=DOFPA+1 ENDDO IF (DOFPA.EQ.0) THEN WRITE(*,*)'0 0 0 SPECIFIES NO MOTION AND IS INVALID' GOTO 111 ENDIF C DO I=1,NSPEC DO J=1,NIONS(I) IF ((RMOVE(I,J).NE.0).AND.(INVERT(I,J).GT.0)) THEN DOF=DOF+DOFPA END IF IF (RMOVE(I,J).NE.0) NMOVE=NMOVE+1 END DO END DO C WRITE(*,*)' WE APPEAR TO HAVE ',DOF,' DEGREES OF FREEDOM' WRITE(*,*)NMOVE,' ATOMS APPEAR TO BE MOVEABLE' C NOFF(1)=0 NIONST=NIONS(1) DO I=2,NSPEC NOFF(I)=NOFF(I-1)+NIONS(I-1) NIONST=NIONST+NIONS(I) END DO C WRITE(*,*)' NIONST=',NIONST C 112 CONTINUE WRITE(*,*)' PLEASE CHOOSE OUTPUT FORMAT' WRITE(*,*)' 1) CETEP' WRITE(*,*)' 2) MSI CASTEP V2' READ(*,*)ICHOICE IF (ICHOICE.LT.1 .OR. ICHOICE.GT.2) GOTO 112 CETEP=.FALSE. IF (ICHOICE.EQ.1) CETEP=.TRUE. C NAT=1 NSP=1 NDOF=1 C IF (CETEP) THEN OPEN (10,FILE='cnstr.dat') WRITE(10,*)DOF ELSE OPEN (10,FILE='c.optim') C No unit cell relaxation WRITE(10,*)0,0,0,0,0,0 C This is the number of non-zero elements in transform marix WRITE(10,*)NMOVE*DOFPA ENDIF C DO I=1,NIONST WRITE(*,*)nsp,nat IF (RMOVE(NSP,NAT).NE.0) THEN IF (INVERT(NSP,NAT).GT.0) THEN C C I.E. (NSP,NAT) IS UNIQUE, OR ITS PAIR HASN'T BEEN CONSIDERED C DO J=1,3 IF (AMOVE(J).NE.0) THEN IF (CETEP) THEN DO K=1,3*NIONST LINE(K)=0 END DO LINE(3*(NOFF(NSP)+INVERT(NSP,NAT)-1)+J)=-1 C INVERT(NSP,NAT)=NAT IFF NO INVERSION PARTNER LINE(3*(NOFF(NSP)+NAT-1)+J)=1 WRITE(10,1111) (LINE(K),K=1,NIONST*3) ELSE WRITE(10,1112) 3*(NOFF(NSP)+NAT-1)+J,NDOF,1.0 IF (INVERT(NSP,NAT).NE.NAT) THEN WRITE(10,1112) & 3*(NOFF(NSP)+INVERT(NSP,NAT)-1)+J,NDOF,-1.0 END IF NDOF=NDOF+1 END IF END IF END DO ENDIF END IF NAT=NAT+1 IF (NAT.GT.NIONS(NSP))THEN NSP=NSP+1 NAT=1 END IF END DO 1111 FORMAT(22I3) 1112 FORMAT(2I7,F6.1) C IF(.NOT.CETEP) THEN do i=1,DOF do j=1,DOF hess(j,i)=0.0 end do end do c do i=1,DOF hess(i,i)=100.0 enddo c write(10,102)((hess(j,i),j=1,i),i=1,DOF) 102 format(6F13.7) ENDIF CLOSE(10) C STOP END