      SUBROUTINE CWRITE
C
C Write selected coordinates on specified file in PDB format.
C for syntax see "HELP"
C
C Author: Axel T. Brunger
C =======================
C
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'coordc.inc'
      INCLUDE 'cnst.inc'
      INCLUDE 'deriv.inc'
C local
      CHARACTER*4 SFROM, FORM
      INTEGER FLAGS, NS
C defaults
      OFILE='OUTPUT'
      FLAGS=ALLHP(INTEG4(NATOM))
      CALL FILL4(HEAP(FLAGS),NATOM,1)
      NS=NATOM
      SFROM='MAIN'
      FORM='CNS'
C
C parsing
      CALL PUSEND('WRITE-COOR>')
      DO WHILE (.NOT.DONE)
      CALL NEXTWD('WRITE-COOR>')
C
      IF (WD(1:4).EQ.'HELP') THEN
C
      CALL CNSHELP('cns-write-coordinate')
C
      ELSE IF (WD(1:4).EQ.'OUTP') THEN
      CALL NEXTFI('OUTPut=',OFILE)
      ELSE IF (WD(1:4).EQ.'FROM') THEN
      CALL NEXTA4('FROM=',SFROM)
      ELSE IF (WD(1:4).EQ.'FORM') THEN
      CALL NEXTA4('FORMat=',FORM)
      ELSE IF (WD(1:4).EQ.'SELE') THEN
      CALL SELCTA(HEAP(FLAGS),NS,X,Y,Z,.TRUE.)
      ELSE
      CALL CHKEND('WRITE-COOR>',DONE)
      END IF
      END DO
      DONE=.FALSE.
C
      IF (SFROM.EQ.'MAIN') THEN
      CALL CWRIT2(OFILE,X,Y,Z,WMAIN,QMAIN,HEAP(FLAGS),NS,FORM)
      ELSE IF (SFROM.EQ.'COMP') THEN
      CALL CWRIT2(OFILE,XCOMP,YCOMP,ZCOMP,WCOMP,QCOMP,HEAP(FLAGS),
     &            NS,FORM)
      ELSE IF (SFROM.EQ.'REFE') THEN
      CALL CWRIT2(OFILE,REFX,REFY,REFZ,KCNSTR,KCNSTR,HEAP(FLAGS),
     &            NS,FORM)
      ELSE IF (SFROM.EQ.'DERI') THEN
      CALL CWRIT2(OFILE,DX,DY,DZ,FBETA,FBETA,HEAP(FLAGS),NS,FORM)
      ELSE
      WRITE(6,'(A)') ' %WRITE-COOR-ERR: unknown FROM parameter:',SFROM
      END IF
C
      CALL FREHP(FLAGS,INTEG4(NATOM))
C
      RETURN
      END
C
      SUBROUTINE CWRIT2(FILE,X,Y,Z,WMAIN,QMAIN,ISLCT,NSLCT,FORM)
C
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'mtf.inc'
      CHARACTER*(*) FILE
      DOUBLE PRECISION X(*), Y(*), Z(*), WMAIN(*), QMAIN(*)
      INTEGER ISLCT(*), NSLCT
      CHARACTER*4 FORM
C local
      INTEGER I, UNIT, L, SEGLEN
      CHARACTER*4 NAME, SEGTMP
      CHARACTER*5 A
      LOGICAL ERROR
C begin
      IF (NSLCT.LE.0) THEN
      WRITE(6,'(A)') ' CWRITE: zero atoms specified. No file created'
      ELSE
      CALL ASSFIL(FILE,UNIT,'WRITE','FORMATTED',ERROR)
      IF (.NOT.ERROR) THEN
      IF(NSLCT.LT.NATOM) WRITE(6,'(A)') ' CWRITE: using atom subset.'
C
C use Brookhaven PDB format (sort of...)
C
CTOM   1223  O   GLY   153     -11.704  -9.200    .489  1.00  0.80
C           XXXX XXXX XXXXy   XXXXXXXXTTTTTTTTUUUUUUUUIIIIIIOOOOOO      !FRODO
C           IIII UUUU XXXXy   OOOOOOOOIIIIIIIIPPPPPPPP      XXXXXX      ;;;;
C
C     ccccc ''''Iyyy O,,,,L   ........>>>>>>>>////////ppppppiiiiii iii  !PDB
C                         ^ insertion character
C                    ^ chain identifier
C           ^ additional character for some atom names (mostly h's)
C
      IF (FORM(1:3).EQ.'CNS'.OR.FORM(1:4).EQ.'PDBO') THEN
      CALL WRTITL(UNIT,1)
      END IF
      DO I=1,NATOM
      IF (ISLCT(I).EQ.1) THEN
C
C right-adjust the resid
      L=4
      CALL TRIMM(RESID(I),L)
      A='     '
      IF (LGE(RESID(I)(L:L),'A')) THEN
C shift it by one unit to the right if an insertion character is
C specified
      A(5-L+1:5)=RESID(I)(1:L)
      ELSE
      A(4-L+1:4)=RESID(I)(1:L)
      END IF
C
CCC      IF (L.EQ.4.OR.LGE(RESID(I)(L:L),'A')) THEN
CCC      A(4-L+1:4)=RESID(I)(1:L)
CCC      ELSE
CCC      A(4-L:3)=RESID(I)(1:L)
CCC      END IF
C
C shift atom names when they exceed 3 characters
      IF (TYPE(I)(4:4).EQ.' ') THEN
      NAME=' '//TYPE(I)(1:3)
      ELSE
      NAME=TYPE(I)
      END IF
C
      IF (FORM.EQ.'PDBA') THEN
C PDB ATOM record
      WRITE(UNIT,
     & '(A,I6,1X,A4,1X,A4,A1,A5,3X,3F8.3,2F6.2,6X,A4)') 'ATOM ',
     & I,NAME,RES(I),IAC(I)(1:1),A,
     & X(I),Y(I),Z(I),QMAIN(I),WMAIN(I),SEGID(I)
      ELSEIF (FORM.EQ.'PDBH') THEN
C PDB HETATM record
      IF (I.GT.99999) THEN
      WRITE(UNIT,
     & '(A,I6,1X,A4,1X,A4,A1,A5,3X,3F8.3,2F6.2,6X,A4)') 'HETAT',
     & I,NAME,RES(I),IAC(I)(1:1),A,
     & X(I),Y(I),Z(I),QMAIN(I),WMAIN(I),SEGID(I)
      ELSE
      WRITE(UNIT,
     & '(A,I5,1X,A4,1X,A4,A1,A5,3X,3F8.3,2F6.2,6X,A4)') 'HETATM',
     & I,NAME,RES(I),IAC(I)(1:1),A,
     & X(I),Y(I),Z(I),QMAIN(I),WMAIN(I),SEGID(I)
      ENDIF
      ELSEIF (FORM.EQ.'PDBO') THEN
C PDB O record
C use segid as chainid if it is one character
      SEGTMP=SEGID(I)
      SEGLEN=4
      CALL TRIMM(SEGTMP,SEGLEN)
      CALL TRIML(SEGTMP,SEGLEN)
      IF (SEGLEN.EQ.1) THEN
      WRITE(UNIT,
     & '(A,I6,1X,A4,1X,A4,A1,A5,3X,3F8.3,2F6.2,6X,A4)') 'ATOM ',
     & I,NAME,RES(I),SEGTMP(1:1),A,
     & X(I),Y(I),Z(I),QMAIN(I),WMAIN(I),SEGID(I)
      ELSE
      WRITE(UNIT,
     & '(A,I6,1X,A4,1X,A4,A1,A5,3X,3F8.3,2F6.2,6X,A4)') 'ATOM ',
     & I,NAME,RES(I),' ',A,
     & X(I),Y(I),Z(I),QMAIN(I),WMAIN(I),SEGID(I)
      END IF
      ELSE
C CNS record, no CHAIN-ID!
      WRITE(UNIT,
     & '(A,I6,1X,A4,1X,A4,1X,A5,3X,3F8.3,2F6.2,6X,A4)') 'ATOM ',
     & I,NAME,RES(I),A,
     & X(I),Y(I),Z(I),QMAIN(I),WMAIN(I),SEGID(I)
      END IF
C
      END IF
      END DO
      IF (FORM(1:3).EQ.'CNS'.OR.FORM(1:4).EQ.'PDBO') THEN
      WRITE(UNIT,'(A)') 'END'
      CALL VCLOSE(UNIT,'KEEP',ERROR)
      END IF
C
      END IF
      END IF
C
      RETURN
      END
C
      SUBROUTINE CREAD(X,Y,Z,WMAIN,QMAIN,FLAGS,QCHAIN)
C
C reads coordinates.
C
C Author: Axel T. Brunger
C
C
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'mtf.inc'
      DOUBLE PRECISION X(*), Y(*), Z(*), WMAIN(*), QMAIN(*)
      INTEGER FLAGS(*)
      LOGICAL QCHAIN
C local
      LOGICAL COND, ECHOLD, QATOM
      INTEGER NMISS, NMULT, I, IATOM, IOFFS
      DOUBLE PRECISION XIN, YIN, ZIN, WIN, QIN
      CHARACTER*4 SID, RID, REN, IUP
      CHARACTER*5 A
      CHARACTER*1 CHAIN
      INTEGER NUMBER, IJ
C parameter
      INTEGER MARK
      DOUBLE PRECISION ONE
      PARAMETER (MARK=-9999, ONE=1.0D0)
C begin
C
C initialize various counters
      NMISS=0
      NMULT=0
      IOFFS=0
      ECHOLD=QECHO
      USED=.TRUE.
      NUMBER=-1
C
      DO WHILE (USED)
      QATOM=.FALSE.
C
      IF (WD(1:4).EQ.'ATOM'.OR.WD(1:4).EQ.'HETA') THEN
C ==================================================
C this section reads the Brookhaven atom coordinates
C ==================================================
      READ(COMLYN,'(12X,A4,1X,A4,A1,A5,3X,3F8.3,F6.2,F6.2,6X,A4)',
     &     ERR=8888) IUP,REN,CHAIN,A,XIN,YIN,ZIN,QIN,WIN,SID
C
C make RESID left-justified and remove trailing blanks
      IJ=5
      CALL TRIML(A,IJ)
      CALL TRIMM(A,IJ)
      IF (IJ.GT.4) THEN
      WRITE(6,'(A,A,A)')
     & ' %CREAD-ERR: residue ID and insertion character ',
     & A,' exceed 4 characters.'
      END IF
      RID=A(1:4)
C
C make IUP left-justified
      IJ=4
      CALL TRIML(IUP,IJ)
C
C make REN left-justified
      IJ=4
      CALL TRIML(REN,IJ)
C
      IF (QCHAIN) THEN
         IF (CHAIN.NE.' ') THEN
            SID=CHAIN
         END IF
      END IF
C
      QATOM=.TRUE.
C
      ELSE IF (WD(1:3).EQ.'TER') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'BREAK') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'CRYST') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:4).EQ.'ORIG') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'SCALE') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:4).EQ.'REMA') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'HEADER') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'COMPND') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'TITLE') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'SOURCE') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'KEYWDS') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'EXPDTA') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'DBREF') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'AUTHOR') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'REVDAT') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:4).EQ.'JRNL') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'SEQADV') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'SEQRES') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'MODRES') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'FTNOTE') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'HET   ') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'HETNAM') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'FORMUL') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'HELIX') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'SHEET') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:4).EQ.'TURN') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'SSBOND') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:4).EQ.'SITE') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:4).EQ.'LINK') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:5).EQ.'MTRIX') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'CISPEP') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'CONECT') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE IF (WD(1:6).EQ.'MASTER') THEN
      QECHO=.TRUE.
      CURSOR=COMLEN
C
      ELSE
      USED=.FALSE.
      END IF
C
      IF (QATOM) THEN
C
C ======================================
C try to find atom in molecular topology
C ======================================
      IOFFS=IOFFS+1
C
C make RID right-justified
      IJ=4
      CALL TRIML(RID,IJ)
C
C get atom number (note the shortcut for coordinate files
C with correct sequence !!)
      COND=SID.EQ.SEGID(IOFFS).AND.RID.EQ.RESID(IOFFS).AND.
     &     REN.EQ.RES(IOFFS).AND.IUP.EQ.TYPE(IOFFS)
      IF (COND) THEN
      IATOM=IOFFS
      ELSE
      IATOM=GETATN(SID,RID,REN,IUP,MARK)
      IOFFS=IATOM
      END IF
      IF (IATOM.EQ.MARK) THEN
      WRITE(6,'(9A)') ' %READC-ERR: atom ',SID,' ',RID,' ',
     &      REN,' ',IUP,' not found in molecular structure'
      ELSE
      IF (FLAGS(IATOM).EQ.1) THEN
      IF (INITIA(IATOM,X,Y,Z)) NMULT=NMULT+1
      X(IATOM)=XIN
      Y(IATOM)=YIN
      Z(IATOM)=ZIN
      QMAIN(IATOM)=QIN
      WMAIN(IATOM)=WIN
      END IF
      END IF
C
      CURSOR=COMLEN
      END IF
C
      IF (USED) THEN
      CALL NEXTWD('COOR>')
      QECHO=.FALSE.
      END IF
C
      END DO
C
C check multiple coordinate assignments
      IF (NMULT.GT.0) THEN
      WRITE(6,'(A,I6,A)') ' %READC-WRN: multiple coordinates for ',
     &     NMULT,' atoms'
      END IF
C
C check still undefined atoms
      IF (IOFFS.GT.0) THEN
      NMISS=0
      DO I=1,NATOM
      IF (.NOT.INITIA(I,X,Y,Z).AND.FLAGS(I).EQ.1) NMISS=NMISS+1
      END DO
      IF (NMISS.GT.0) THEN
      WRITE(6,'(A,I6,A)') ' %READC-WRN: still ',NMISS,
     &              ' missing coordinates (in selected subset)'
      END IF
      END IF
C
C-ERROR-LABEL
      GOTO 7777
8888  ERROR=.TRUE.
      WRITE(6,'(A)') ' %COOR-ERR: ERROR during reading coordinates'
7777  CONTINUE
C
      QECHO=ECHOLD
      RETURN
      END
C
