text:
      SUBROUTINE SELCTA(FLAGS,NSELCT,X,Y,Z,QCOOR)
C
C Recursive atom selection routine.
C
C For syntax see main parsing loop "HELP"
C
C Author: Axel T. Brunger
C =======================
C
C Parameter MAXDEP specifies the maximum depth of an expression
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'mtf.inc'
      INTEGER FLAGS(*)
      INTEGER NSELCT
      DOUBLE PRECISION X(*), Y(*), Z(*)
      LOGICAL QCOOR
C local
      INTEGER RETADR, QSTACK
      INTEGER RETMAX, MAXDEP, STKSIZ
      PARAMETER (RETMAX=100, MAXDEP=10)
C begin
      STKSIZ=MAX(MAXDEP*NATOM,MAXDEP)
      RETADR=ALLHP(INTEG4(RETMAX))
      QSTACK=ALLHP(ILOGIC(STKSIZ))
      CALL SELCT2(RETMAX,HEAP(RETADR),FLAGS,NSELCT,QCOOR,X,Y,Z,
     &            STKSIZ,HEAP(QSTACK))
      CALL FREHP(QSTACK,ILOGIC(STKSIZ))
      CALL FREHP(RETADR,INTEG4(RETMAX))
      END
C
      SUBROUTINE SELCT2(RETMAX,RETADR,FLAGS,NSELCT,QCOOR,
     &           XC,YC,ZC,STKSIZ,QSTACK)
C
C Author: Axel T. Brunger
C =======================
C
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'vector.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'symbol.inc'
      INCLUDE 'timer.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'coordc.inc'
      INCLUDE 'deriv.inc'
      INCLUDE 'cnst.inc'
      INCLUDE 'consta.inc'
      INCLUDE 'xcrystal.inc'
      INCLUDE 'xtarget.inc'
      INCLUDE 'ncs.inc'
      INTEGER RETMAX, RETADR(*), FLAGS(*), NSELCT
      LOGICAL QCOOR
      DOUBLE PRECISION XC(*), YC(*), ZC(*)
      INTEGER STKSIZ
      LOGICAL QSTACK(*)
C local
      INTEGER QSTKLN, QSTART, QINDX1, QINDX2, QLIFT
      INTEGER ISTK
      LOGICAL ELEMNT, OK, ERR, QKNOWN, XRERR
      LOGICAL QRANGE, QABS, QWCARD, QWCARD1, QWCARD2, QWCARD3, EXIST
      LOGICAL CLOOP
      INTEGER IAT, IATS1
      INTEGER IAT2, NUMB, NUMB1, NUMB2
      INTEGER I, INDX1, IGRP, IBT, JBT
      INTEGER J, IDLEN, K, IATS2, IIAT
      INTEGER LEVEL, ISYM, INSYM
      DOUBLE PRECISION RMAX, RMAX2, DIST2, DXL, DYL, DZL, R(3)
      DOUBLE PRECISION XREF, YREF, ZREF
C MODIFICATION: selection token lengths doubled for wildcards with
C backslash escapes, and for input of numeric resSeq > 9999.
      CHARACTER*8 SEG1,  RES1,  RESN1,  RESID1,  SEG2,  RESN2
      INTEGER LSEG1, LRESN1, LRESD1, LSEG2, LRESN2
      CHARACTER*8 RESID2, TYPE1, TYPE2
      INTEGER LRESD2,LTYPE1,LTYPE2
      INTEGER LRESD
      CHARACTER*2 ALT1, ALT2
      INTEGER LALT1, LALT2
      CHARACTER*1 ALPH, ALPH1, ALPH2
      CHARACTER*1 SBRA, SKET
      PARAMETER (SBRA='(',SKET=')')
      DOUBLE PRECISION DPVAL
      DOUBLE COMPLEX DCVAL
      DOUBLE PRECISION XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX
      DOUBLE PRECISION XRRF, YRRF, ZRRF, XRRFF, YRRFF, ZRRFF
C pointer
      INTEGER XF, YF, ZF, XFF, YFF, ZFF, XFFF, YFFF, ZFFF
C parameter
      DOUBLE PRECISION ZERO, ONE, BIG
      PARAMETER (ZERO=0.0D0, ONE=1.0D0, BIG=10000.0D0)
C
C begin
C
C initialize working stack
      QSTART=1
      QSTKLN=NATOM
      QLIFT=0
C
      LEVEL=1
C
      CALL SENEXT
      IF (WD(1:4).EQ.'=   ') CALL SENEXT
      IF (WD(1:2).EQ.'? ') THEN
C
C-begin-show-the-current-set
C
C This procedure prints a listing of the atoms that are currently
C set.
      WRITE(6,7001)
7001  FORMAT(' The following atoms are currently set:')
C
      SEG1='    '
      RES1='    '
      RESN1='    '
      J=0
      DO IAT=1,NATOM
      IF (SEG1.NE.SEGID(IAT).OR.RES1.NE.RESID(IAT).OR.RESN1.NE.RES(IAT))
     1 THEN
      SEG1=SEGID(IAT)
      RES1=RESID(IAT)
      RESN1=RES(IAT)
      IF (J.GT.0) WRITE(6,'(2A)') ' ',WD(1:WDLEN)
      WDLEN=0
      CALL ADDST(WD,WDMAX,WDLEN,SEGID(IAT),4)
      CALL ADDST(WD,WDMAX,WDLEN,' ',1)
      CALL ADDST(WD,WDMAX,WDLEN,RESID(IAT),4)
      CALL ADDST(WD,WDMAX,WDLEN,' ',1)
      CALL ADDST(WD,WDMAX,WDLEN,RES(IAT),4)
      IDLEN = WDLEN
      J=0
      END IF
      IF (FLAGS(IAT).EQ.1) THEN
      IF (WDLEN+1+4 .GT. 78) THEN
      WRITE(6,'(2A)') ' ',WD(1:WDLEN)
      J = 0
      WDLEN = 0
      DO K=1,IDLEN
      CALL ADDST(WD,WDMAX,WDLEN,' ',1)
      END DO
      END IF
      CALL ADDST(WD,WDMAX,WDLEN,' ',1)
      CALL ADDST(WD,WDMAX,WDLEN,TYPE(IAT),4)
      J=J+1
      END IF
      END DO
      IF (J.GT.0) WRITE(6,'(2A)') ' ',WD(1:WDLEN)
C-end-show-the-current-set
C
C jump to exit
      GOTO 99999
C
      END IF
C
C check that expression begins with a "("
      IF (.NOT.WD(1:1).EQ.SBRA) THEN
      CALL DSPERR('SELRPN','selection has to begin with a "("')
      CALL WRNDIE(-5,'SELRPN','parsing error')
C jump to exit
      GOTO 99999
      END IF
C
      CALL SENEXT
C ** invoke procedure expression **
      RETADR(LEVEL)=1
      GOTO 1110
1111  CONTINUE
C ** return label                **
C
C
      IF (.NOT.WD(1:1).EQ.SKET) THEN
      CALL DSPERR('SELRPN','selection has to end with a ")"')
      CALL SAVEWD
      END IF
      DONE=.FALSE.
C
C ** return      **
      GOTO 9999
C ** return      **
C
C
C---- BEGIN PROCEDURE EXPRESSION --------------------------------------
1110  CONTINUE
      CALL SEDOWN(LEVEL,RETMAX)
C
C ** invoke procedure term    **
      RETADR(LEVEL)=3
      GOTO 2220
2223  CONTINUE
C ** return label             **
C
2233  IF (WD(1:2).NE.'OR') GOTO 1122
      CALL SENEXT
C
C ** invoke procedure term    **
      RETADR(LEVEL)=4
      GOTO 2220
2224  CONTINUE
C ** return label             **
C
C-begin-process-or-operator
      QINDX1=QSTART
      QINDX2=QSTART-QSTKLN
      DO IAT=1,NATOM
      QSTACK(QINDX2)=QSTACK(QINDX1).OR.QSTACK(QINDX2)
      QINDX1=QINDX1+1
      QINDX2=QINDX2+1
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
C-end-process-or-operator
C
      GOTO 2233
1122  CONTINUE
C
C
      CALL SEUP(LEVEL,RETMAX)
C
C     return to address:
      IF (RETADR(LEVEL).EQ.1) GOTO 1111
      IF (RETADR(LEVEL).EQ.2) GOTO 1112
      WRITE(6,'(A)')
     & ' %SELRPN-ERR: Return address unknown. Fatal coding error'
      CALL DIE
C---- END PROCEDURE EXPRESSION ---------------------------------------
C
C
C
C---- BEGIN PROCEDURE TERM -------------------------------------------
2220  CONTINUE
      CALL SEDOWN(LEVEL,RETMAX)
C
C ** invoke procedure factor   **
      RETADR(LEVEL)=5
      GOTO 3330
3335  CONTINUE
C ** return label              **
C
4455  IF (WD(1:3).NE.'AND') GOTO 5566
      CALL SENEXT
C
C ** invoke procedure factor   **
      RETADR(LEVEL)=6
      GOTO 3330
3336  CONTINUE
C ** return label              **
C
C-begin-process-and-operator
      QINDX1=QSTART
      QINDX2=QSTART-QSTKLN
      DO IAT=1,NATOM
      QSTACK(QINDX2)=QSTACK(QINDX1).AND.QSTACK(QINDX2)
      QINDX1=QINDX1+1
      QINDX2=QINDX2+1
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
C-end-process-and-operator
C
      GOTO 4455
5566  CONTINUE
C
      CALL SEUP(LEVEL,RETMAX)
C
C     return to address:
      IF (RETADR(LEVEL).EQ.3) GOTO 2223
      IF (RETADR(LEVEL).EQ.4) GOTO 2224
      WRITE(6,'(A)')
     & ' %SELRPN-ERR: Return address unknown. Fatal coding error.'
      CALL DIE
C---- END PROCEDURE TERM ----------------------------------------------
C
C
C
C
C
C---- BEGIN PROCEDURE FACTOR ------------------------------------------
3330  CONTINUE
      CALL SEDOWN(LEVEL,RETMAX)
C=NOT===================================================================
      IF (WD(1:3).NE.'NOT') GOTO 10
C
      CALL SENEXT
C
C ** invoke procedure factor  **
      RETADR(LEVEL)=7
      GOTO 3330
3337  CONTINUE
C ** return label             **
C
C-begin-process-not-operator
      QINDX1=QSTART
      DO IAT=1,NATOM
      QSTACK(QINDX1)=.NOT.QSTACK(QINDX1)
      QINDX1=QINDX1+1
      END DO
C-end-process-not-operator
C
      GOTO 77777
C=BYRE===================================================================
10    IF (WD(1:4).NE.'BYRE') GOTO 20
C
      CALL SENEXT
C
C ** invoke procedure factor  **
      RETADR(LEVEL)=8
      GOTO 3330
3338  CONTINUE
C ** return label             **
C
C-begin-process-byres-operator
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1)=QSTACK(INDX1-QSTKLN)
      END DO
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1-QSTKLN)=.FALSE.
      END DO
      DO IAT2=1,NATOM
      IF (QSTACK(QSTART+IAT2-1)) THEN
      SEG1=SEGID(IAT2)
      RES1=RESID(IAT2)
      RESN1=RES(IAT2)
      DO IAT=1,NATOM
      IF (SEGID(IAT).EQ.SEG1) THEN
      IF (RESID(IAT).EQ.RES1.AND.RES(IAT).EQ.RESN1) THEN
      QSTACK(QSTART-QSTKLN+IAT-1)=.TRUE.
      END IF
      END IF
      END DO
      END IF
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
C-end-process-byres-operator
C
      GOTO 77777
C=BYGR==============================================================
20    IF (WD(1:4).NE.'BYGR') GOTO 30
C
      CALL SENEXT
C
C ** invoke procedure factor  **
      RETADR(LEVEL)=9
      GOTO 3330
3339  CONTINUE
C ** return label             **
C
C-begin-process-bygroup-operator
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1)=QSTACK(INDX1-QSTKLN)
      END DO
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1-QSTKLN)=.FALSE.
      END DO
      DO IGRP=1,NGRP
      IAT2=IGPBS(IGRP)
      ELEMNT=.FALSE.
      CLOOP = .TRUE.
      DO WHILE (CLOOP)
      IAT2=IAT2+1
      IF (QSTACK(QSTART+IAT2-1)) THEN
      DO IAT=IGPBS(IGRP)+1,IGPBS(IGRP+1)
      QSTACK(QSTART-QSTKLN+IAT-1)=.TRUE.
      END DO
      ELEMNT=.TRUE.
      END IF
      IF (ELEMNT.OR.IAT2.GE.IGPBS(IGRP+1)) CLOOP = .FALSE.
      END DO
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
C-end-process-bygroup-operator
C
      GOTO 77777
C=BOUN=========================================================
30    IF (WD(1:4).NE.'BOND') GOTO 40
C
      CALL SENEXT
C
C ** invoke procedure factor  **
      RETADR(LEVEL)=10
      GOTO 3330
3340  CONTINUE
C ** return label             **
C
C BONDedto operator added by Michael Nilges, EMBL
C-begin-process-bondedto-operator
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1)=QSTACK(INDX1-QSTKLN)
      END DO
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1-QSTKLN)=.FALSE.
      END DO
      DO I=1,NBOND
      IBT=IB(I)
      JBT=JB(I)
      IF (IBT.GT.0 .AND. JBT.GT.0) THEN
      IF(QSTACK(QSTART+JBT-1)) THEN
      QSTACK(QSTART-QSTKLN+IBT-1)=.TRUE.
      END IF
      IF(QSTACK(QSTART+IBT-1)) THEN
      QSTACK(QSTART-QSTKLN+JBT-1)=.TRUE.
      END IF
      END IF
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
C-end-process-boundto-operator
C
      GOTO 77777
C=(============================================================
40    IF (WD(1:1).NE.SBRA) GOTO 50
C
      CALL SENEXT
C
C ** invoke procedure expression **
      RETADR(LEVEL)=2
      GOTO 1110
1112  CONTINUE
C ** return label                **
      IF (.NOT.WD(1:1).EQ.SKET) THEN
      WRITE(6,'(A)') ' %SELRPN-ERR: unbalanced parentheses.'
      CALL WRNDIE(-5,'SELRPN','parsing error')
C exit selection subroutine
      GOTO 99999
      END IF
      CALL SENEXT
      GOTO 77777
C==============================================================
C==========
C ALL token
C==========
50    IF (WD(1:3).EQ.'ALL') THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      ELEMNT=.TRUE.
      QSTACK(QSTART+IAT-1)=ELEMNT
      END DO
      CALL SENEXT
C==========
C NONE token
C==========
      ELSE IF (WD(1:4).EQ.'NONE') THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      ELEMNT=.FALSE.
      QSTACK(QSTART+IAT-1)=ELEMNT
      END DO
      CALL SENEXT
C
C===========
C PREV token
C===========
      ELSE IF (WD(1:4).EQ.'PREV') THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      ELEMNT=(FLAGS(IAT).EQ.1)
      QSTACK(QSTART+IAT-1)=ELEMNT
      END DO
      CALL SENEXT
C
C=============
C RECALL token
C=============
      ELSE IF (WD(1:4).EQ.'RECA'.OR.WD(1:4).EQ.'STOR') THEN
C fill default (FALSE)
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
C
      IF (WD(1:7).EQ.'RECALL1'.OR.WD(1:6).EQ.'STORE1') THEN
      J=1
      ELSE IF (WD(1:7).EQ.'RECALL2'.OR.WD(1:6).EQ.'STORE2') THEN
      J=2
      ELSE IF (WD(1:7).EQ.'RECALL3'.OR.WD(1:6).EQ.'STORE3') THEN
      J=3
      ELSE IF (WD(1:7).EQ.'RECALL4'.OR.WD(1:6).EQ.'STORE4') THEN
      J=4
      ELSE IF (WD(1:7).EQ.'RECALL5'.OR.WD(1:6).EQ.'STORE5') THEN
      J=5
      ELSE IF (WD(1:7).EQ.'RECALL6'.OR.WD(1:6).EQ.'STORE6') THEN
      J=6
      ELSE IF (WD(1:7).EQ.'RECALL7'.OR.WD(1:6).EQ.'STORE7') THEN
      J=7
      ELSE IF (WD(1:7).EQ.'RECALL8'.OR.WD(1:6).EQ.'STORE8') THEN
      J=8
      ELSE IF (WD(1:7).EQ.'RECALL9'.OR.WD(1:6).EQ.'STORE9') THEN
      J=9
      ELSE
      CALL NEXTI('RECALL=',J)
      END IF
      IF (.NOT. (J.GE.1.AND.J.LE.MAXSTO)) THEN
      CALL WRNDIE(-5,'SELRPN','Store number out of range')
      ELSE IF (PTRSTO(J).EQ.0) THEN
      WRITE(6,'(A)') ' SELRPN: Store empty. No atoms selected.'
      ELSE IF (LENSTO(J).NE.NATOM) THEN
      CALL WRNDIE(-5,'SELRPN','dim. mismatch')
      ELSE
      CALL SELRCL(QSTACK(QSTART),HEAP(PTRSTO(J)),NATOM)
      END IF
      CALL SENEXT
C
C=========
C ID token
C=========
      ELSE IF (WD(1:4).EQ.'BYNU'.OR.WD(1:2).EQ.'ID') THEN
      CALL SENEXT
      IATS1=DECODI(WD,WDLEN,OK)
      ERR=.FALSE.
      IF (.NOT.OK.OR.IATS1.LE.0.OR.IATS1.GT.NATOM) THEN
      WRITE(6,'(2A)') ' %SELRPN-BYNU-ERR: no valid number ',WD(1:WDLEN)
      CALL WRNDIE(-5,'SELRPN','no valid number')
      ERR=.TRUE.
      END IF
      IATS2=IATS1
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL SENEXT
      IATS2=DECODI(WD,WDLEN,OK)
      IF (.NOT.OK.OR.IATS2.LE.0.OR.IATS2.GT.NATOM) THEN
      WRITE(6,'(2A)') ' %SELRPN-BYNU-ERR: no valid number ',WD(1:WDLEN)
      CALL WRNDIE(-5,'SELRPN','no valid number')
      ERR=.TRUE.
      END IF
      CALL SENEXT
      END IF
      IF (.NOT.ERR) THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
      DO IAT=IATS1,IATS2
      QSTACK(QSTART+IAT-1)=.TRUE.
      END DO
      END IF
C
C============
C SEGID token
C============
      ELSE IF (WD(1:4).EQ.'SEGI') THEN
      CALL NEXTST('segment-id=',WD)
C
      QWCARD=.NOT.QQUOT
C
      CALL COPYST(SEG1,LEN(SEG1),LSEG1,WD,WDLEN)
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL NEXTST('segment-id=',WD)
      CALL COPYST(SEG2,LEN(SEG2),LSEG2,WD,WDLEN)
      CALL SENEXT
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=
     &    (SEG1(1:LSEG1).LE.SEGID(IAT).AND.SEGID(IAT).LE.SEG2(1:LSEG2))
      END DO
      ELSE
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD) THEN
      CALL EQSTWC(SEGID,4,SEG1,LSEG1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=SEG1(1:LSEG1).EQ.SEGID(IAT)
      END DO
      END IF
C
      END IF
C
C============
C ALTID token
C============
      ELSE IF (WD(1:4).EQ.'ALTI') THEN
      CALL NEXTST('alternate-id=',WD)
C
      QWCARD=.NOT.QQUOT
C
      CALL COPYST(ALT1,LEN(ALT1),LALT1,WD,WDLEN)
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL NEXTST('alternate-id=',WD)
      CALL COPYST(ALT2,LEN(ALT1),LALT2,WD,WDLEN)
      CALL SENEXT
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=
     &    (ALT1(1:LALT1).LE.ALTID(IAT).AND.ALTID(IAT).LE.ALT2(1:LALT2))
      END DO
      ELSE
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD) THEN
      CALL EQSTWC(ALTID,1,ALT1,LALT1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=ALT1(1:LALT1).EQ.ALTID(IAT)
      END DO
      END IF
C
      END IF
C
C===========
C RESN token
C===========
      ELSE IF (WD(1:4).EQ.'RESN') THEN
      CALL NEXTST('residue-name=',WD)
C
      QWCARD=.NOT.QQUOT
C
      CALL COPYST(RESN1,LEN(RESN1),LRESN1,WD,WDLEN)
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL NEXTST('residue-name=',WD)
      CALL COPYST(RESN2,LEN(RESN2),LRESN2,WD,WDLEN)
      CALL SENEXT
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=
     &   (RESN1(1:LRESN1).LE.RES(IAT).AND.RES(IAT).LE.RESN2(1:LRESN2))
      END DO
      ELSE
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD) THEN
      CALL EQSTWC(RES,4,RESN1,LRESN1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=RESN1(1:LRESN1).EQ.RES(IAT)
      END DO
      END IF
C
      END IF
C
C============
C RESID token
C============
      ELSE IF (WD(1:4).EQ.'RESI') THEN
      CALL NEXTST('residue-id=',WD)
C
      QWCARD=.NOT.QQUOT
C
      CALL COPYST(RESID1,LEN(RESID1),LRESD1,WD,WDLEN)
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL SPLITI(RESID1,NUMB1,ALPH1,LRESD1,OK)
C
      IF (.NOT.OK) THEN
      WRITE(6,'(2A)') ' %SELRPN-ERR: invalid RESID ',RESID1
      CALL WRNDIE(-5,'SELRPN','invalid residue specification')
      ERR=.TRUE.
      END IF
C
      CALL NEXTST('residue-id=',WD)
      CALL COPYST(RESID2,LEN(RESID2),LRESD2,WD,WDLEN)
      CALL SPLITI(RESID2,NUMB2,ALPH2,LRESD2,OK)
C
      IF (.NOT.OK) THEN
      WRITE(6,'(2A)') ' %SELRPN-ERR: invalid RESID ',RESID2
      CALL WRNDIE(-5,'SELRPN','invalid residue specification')
      ERR=.TRUE.
      END IF
C
      CALL SENEXT
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      LRESD=4
      CALL SPLITI(RESID(IAT),NUMB,ALPH,LRESD,OK)
C
      IF (.NOT.OK) THEN
      WRITE(6,'(3A,I9)')
     & ' %SELRPN-ERR: invalid RESId ',RESID(IAT),' for atom ',IAT
      CALL WRNDIE(-5,'SELRPN','invalid RESID.')
      ERR=.TRUE.
      END IF
C
      QSTACK(QSTART+IAT-1)=
     &      (NUMB1.LT.NUMB.OR.
     &      (NUMB1.EQ.NUMB.AND.ALPH1.LE.ALPH))
     &       .AND.
     &      (NUMB.LT.NUMB2.OR.
     &      (NUMB.EQ.NUMB2.AND.ALPH.LE.ALPH2))
      END DO
      ELSE
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD) THEN
      CALL EQSTWC(RESID,4,RESID1,LRESD1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=RESID1(1:LRESD1).EQ.RESID(IAT)
      END DO
      END IF
C
      END IF
C
C===========
C NAME token
C===========
      ELSE IF (WD(1:4).EQ.'TYPE'.OR.WD(1:4).EQ.'NAME') THEN
      CALL NEXTST('atom-name=',WD)
C
      QWCARD=.NOT.QQUOT
C
      CALL COPYST(TYPE1,LEN(TYPE1),LTYPE1,WD,WDLEN)
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL NEXTST('atom-name=',WD)
      CALL COPYST(TYPE2,LEN(TYPE2),LTYPE2,WD,WDLEN)
      CALL SENEXT
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=
     &  TYPE1(1:LTYPE1).LE.TYPE(IAT).AND.TYPE(IAT).LE.TYPE2(1:LTYPE2)
      END DO
      ELSE
C
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD) THEN
      CALL EQSTWC(TYPE,4,TYPE1,LTYPE1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=TYPE1(1:LTYPE1).EQ.TYPE(IAT)
      END DO
      END IF
C
      END IF
C
C===========
C CHEM token
C===========
      ELSE IF (WD(1:4).EQ.'CHEM') THEN
      CALL NEXTST('chemical-type=',WD)
C
      QWCARD=.NOT.QQUOT
C
      CALL COPYST(TYPE1,LEN(TYPE1),LTYPE1,WD,WDLEN)
      CALL SENEXT
      QRANGE=(WD(1:1).EQ.':')
      IF (QRANGE) THEN
      CALL NEXTST('chemical-type=',WD)
      CALL COPYST(TYPE2,LEN(TYPE2),LTYPE2,WD,WDLEN)
      CALL SENEXT
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=
     &   TYPE1(1:LTYPE1).LE.IAC(IAT).AND.IAC(IAT).LE.TYPE2(1:LTYPE2)
      END DO
      ELSE
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD) THEN
      CALL EQSTWC(IAC,4,TYPE1,LTYPE1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=TYPE1(1:LTYPE1).EQ.IAC(IAT)
      END DO
      END IF
C
      END IF
C
C===========
C ATOM token
C===========
      ELSE IF (WD(1:4).EQ.'ATOM') THEN
      CALL NEXTST('segment-id=',WD)
C
      QWCARD1=.NOT.QQUOT
C
      CALL COPYST(SEG1,LEN(SEG1),LSEG1,WD,WDLEN)
      CALL NEXTST('residue-id=',WD)
C
      QWCARD2=.NOT.QQUOT
C
      CALL COPYST(RES1,LEN(RES1),LRESD1,WD,WDLEN)
      CALL NEXTST('atom-name=',WD)
C
      QWCARD3=.NOT.QQUOT
C
      CALL COPYST(WDT,WDTMAX,WDTLEN,WD,WDLEN)
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD1) THEN
      CALL EQSTWC(SEGID,4,SEG1,LSEG1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=SEG1(1:LSEG1).EQ.SEGID(IAT)
      END DO
      END IF
C
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD2) THEN
      CALL EQSTWC(RESID,4,RES1,LRESD1,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=RES1(1:LRESD1).EQ.RESID(IAT)
      END DO
      END IF
C
      QINDX1=QSTART
      QINDX2=QSTART-QSTKLN
      DO IAT=1,NATOM
      QSTACK(QINDX2)=QSTACK(QINDX1).AND.QSTACK(QINDX2)
      QINDX1=QINDX1+1
      QINDX2=QINDX2+1
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
      IF (QWCARD3) THEN
      CALL EQSTWC(TYPE,4,WDT,WDTLEN,QSTART,NATOM,QSTACK)
      ELSE
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=WDT(1:WDTLEN).EQ.TYPE(IAT)
      END DO
      END IF
C
      QINDX1=QSTART
      QINDX2=QSTART-QSTKLN
      DO IAT=1,NATOM
      QSTACK(QINDX2)=QSTACK(QINDX1).AND.QSTACK(QINDX2)
      QINDX1=QINDX1+1
      QINDX2=QINDX2+1
      END DO
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
      CALL SENEXT
C
C===========
C POIN token
C===========
      ELSE IF (WD(1:4).EQ.'POIN') THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      CALL SELTVF('POINt-vector=',R)
      XREF=R(1)
      YREF=R(2)
      ZREF=R(3)
      CALL SENEXT
      IF (WD(1:3).EQ.'CUT') THEN
      CALL NEXTF('CUT=',RMAX)
      CALL SENEXT
      ELSE
      RMAX=8.0D0
      END IF
      IF (.NOT.QCOOR) THEN
      CALL WRNDIE(-5,'SELRPN','no coordinates present')
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
      ELSE
      RMAX2=RMAX*RMAX
      DO IAT=1,NATOM
      DXL=XREF-XC(IAT)
      DYL=YREF-YC(IAT)
      DZL=ZREF-ZC(IAT)
      DIST2=DXL*DXL+DYL*DYL+DZL*DZL
      ELEMNT=(DIST2.LE.RMAX2)
      QSTACK(QSTART+IAT-1)=ELEMNT
      END DO
      END IF
C
C===========
C ATTR token
C===========
      ELSE IF (WD(1:4).EQ.'ATTR') THEN
      CALL NEXTQL('ATTRibute=')
      IF (WD(1:3).EQ.'ABS') THEN
      QABS=.TRUE.
      CALL NEXTQL('ATTRibute=')
      ELSE
      QABS=.FALSE.
      END IF
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      IF (WD(1:4).EQ.'X   ') THEN
      CALL SELATT(X,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'Y   ') THEN
      CALL SELATT(Y,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'Z   ') THEN
      CALL SELATT(Z,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'WMAI') THEN
      CALL SELATT(WMAIN,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'QMAI') THEN
      CALL SELATT(QMAIN,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'B   ') THEN
      CALL SELATT(WMAIN,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'Q   ') THEN
      CALL SELATT(QMAIN,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'XCOM') THEN
      CALL SELATT(XCOMP,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'YCOM') THEN
      CALL SELATT(YCOMP,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'ZCOM') THEN
      CALL SELATT(ZCOMP,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'WCOM') THEN
      CALL SELATT(WCOMP,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'BCOM') THEN
      CALL SELATT(WCOMP,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'QCOM') THEN
      CALL SELATT(QCOMP,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'REFX') THEN
      CALL SELATT(REFX,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'REFY') THEN
      CALL SELATT(REFY,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'REFZ') THEN
      CALL SELATT(REFZ,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'RMSD') THEN
      CALL SELATT(RMSD,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'HARM') THEN
      CALL SELATT(KCNSTR,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'ZHAR') THEN
      CALL SELATT(KZCNSTR,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'MASS') THEN
      CALL SELATT(AMASS,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'CHAR') THEN
      CALL SELATT(CG,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE1') THEN
      CALL SELATT(HEAP(PTRSTO(1)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE2') THEN
      CALL SELATT(HEAP(PTRSTO(2)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE3') THEN
      CALL SELATT(HEAP(PTRSTO(3)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE4') THEN
      CALL SELATT(HEAP(PTRSTO(4)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE5') THEN
      CALL SELATT(HEAP(PTRSTO(5)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE6') THEN
      CALL SELATT(HEAP(PTRSTO(6)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE7') THEN
      CALL SELATT(HEAP(PTRSTO(7)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE8') THEN
      CALL SELATT(HEAP(PTRSTO(8)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:6).EQ.'STORE9') THEN
      CALL SELATT(HEAP(PTRSTO(9)),QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'DX  ') THEN
      CALL SELATT(DX,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'DY  ') THEN
      CALL SELATT(DY,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'DZ  ') THEN
      CALL SELATT(DZ,QSTART,NATOM,QSTACK,QABS)
      ELSE IF (WD(1:4).EQ.'VX  ') THEN
      DO I=1,NATOM
      XV(I)=XV(I)/TIMFAC
      END DO
      CALL SELATT(XV,QSTART,NATOM,QSTACK,QABS)
      DO I=1,NATOM
      XV(I)=XV(I)*TIMFAC
      END DO
      ELSE IF (WD(1:4).EQ.'VY  ') THEN
      DO I=1,NATOM
      YV(I)=YV(I)/TIMFAC
      END DO
      CALL SELATT(YV,QSTART,NATOM,QSTACK,QABS)
      DO I=1,NATOM
      YV(I)=YV(I)*TIMFAC
      END DO
      ELSE IF (WD(1:4).EQ.'VZ  ') THEN
      DO I=1,NATOM
      ZV(I)=ZV(I)/TIMFAC
      END DO
      CALL SELATT(ZV,QSTART,NATOM,QSTACK,QABS)
      DO I=1,NATOM
      ZV(I)=ZV(I)*TIMFAC
      END DO
      ELSE IF (WD(1:4).EQ.'FBET') THEN
      CALL SELATT(FBETA,QSTART,NATOM,QSTACK,QABS)
C
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_A1') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('A1',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_A2') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('A2',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_A3') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('A3',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_A4') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('A4',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_B1') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('B1',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_B2') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('B2',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_B3') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('B3',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_B4') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('B4',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_C') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('C',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_FP') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('FP',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE IF (WD(1:WDLEN).EQ.'SCATTER_FDP') THEN
      XF=ALLHP(IREAL8(NATOM))
      CALL XSCAFIL('FDP',NATOM,HEAP(XF),XRERR)
      CALL SELATT(HEAP(XF),QSTART,NATOM,QSTACK,QABS)
      CALL FREHP(XF,IREAL8(NATOM))
C
      ELSE
      CALL DSPERR('SELRPN',' unknown attribute')
      END IF
      CALL SENEXT
C
C============
C KNOWn token
C============
      ELSE IF (WD(1:4).EQ.'INIT'.OR.WD(1:4).EQ.'KNOW') THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      ELEMNT=INITIA(IAT,X,Y,Z)
      QSTACK(QSTART+IAT-1)=ELEMNT
      END DO
      CALL SENEXT
C
C===========
C HYDR token
C===========
      ELSE IF (WD(1:4).EQ.'HYDR') THEN
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO IAT=1,NATOM
      ELEMNT=HYDROG(IAT)
      QSTACK(QSTART+IAT-1)=ELEMNT
      END DO
      CALL SENEXT
C
C===========
C TAG token
C===========
      ELSE IF (WD(1:4).EQ.'TAG ') THEN
C
C selects exactly one atom from each "entity".  An
C entity is defined as a set of atoms with equal segids,
C resids, and residue names.
C
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      SEG1='    '
      RES1='    '
      RESN1='    '
      DO IAT=1,NATOM
      IF (SEG1.NE.SEGID(IAT).OR.RES1.NE.RESID(IAT).OR.
     &      RESN1.NE.RES(IAT)) THEN
      SEG1=SEGID(IAT)
      RES1=RESID(IAT)
      RESN1=RES(IAT)
C
C check if it has been previously already defined
      EXIST=.FALSE.
      DO IIAT=1,IAT-1
      IF (QSTACK(QSTART+IIAT-1)) THEN
      IF (SEG1.EQ.SEGID(IIAT).AND.RES1.EQ.RESID(IIAT).AND.
     &      RESN1.EQ.RES(IIAT)) THEN
      EXIST=.TRUE.
      END IF
      END IF
      END DO
C
C
      IF (.NOT.EXIST) THEN
      QSTACK(QSTART+IAT-1)=.TRUE.
      ELSE
      QSTACK(QSTART+IAT-1)=.FALSE.
      END IF
      ELSE
      QSTACK(QSTART+IAT-1)=.FALSE.
      END IF
      END DO
      CALL SENEXT
C=====================================
      ELSE IF (WD(1:4).EQ.'FBOX') THEN
      CALL NEXTF('xmin=',XMIN)
      CALL NEXTF('xmax=',XMAX)
      CALL NEXTF('ymin=',YMIN)
      CALL NEXTF('ymax=',YMAX)
      CALL NEXTF('zmin=',ZMIN)
      CALL NEXTF('zmax=',ZMAX)
      XMIN=MAX(ZERO,XMIN)
      XMAX=MIN(ONE,XMAX)
      YMIN=MAX(ZERO,YMIN)
      YMAX=MIN(ONE,YMAX)
      ZMIN=MAX(ZERO,ZMIN)
      ZMAX=MIN(ONE,ZMAX)
      CALL SENEXT
C
      IF (.NOT.QCOOR) THEN
      CALL WRNDIE(-5,'SELRPN','no coordinates present')
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
      ELSE
      DO IAT=1,NATOM
C
      QSTACK(QSTART+IAT-1)=.FALSE.
C compute fractional coordinates
      CALL XPSYMO(1,XC(IAT),YC(IAT),ZC(IAT),1,XRRFF,YRRFF,ZRRFF)
C
C apply mod operation
      XRRFF=MOD(XRRFF+BIG,ONE)
      YRRFF=MOD(YRRFF+BIG,ONE)
      ZRRFF=MOD(ZRRFF+BIG,ONE)
      IF          (XRRFF.GE.XMIN.AND.XRRFF.LE.XMAX
     &        .AND.YRRFF.GE.YMIN.AND.YRRFF.LE.YMAX
     &        .AND.ZRRFF.GE.ZMIN.AND.ZRRFF.LE.ZMAX) THEN
      QSTACK(QSTART+IAT-1)=.TRUE.
      END IF
      END DO
      END IF
C=====================================
      ELSE IF (WD(1:5).EQ.'SFBOX') THEN
      CALL NEXTF('xmin=',XMIN)
      CALL NEXTF('xmax=',XMAX)
      CALL NEXTF('ymin=',YMIN)
      CALL NEXTF('ymax=',YMAX)
      CALL NEXTF('zmin=',ZMIN)
      CALL NEXTF('zmax=',ZMAX)
      XMIN=MAX(ZERO,XMIN)
      XMAX=MIN(ONE,XMAX)
      YMIN=MAX(ZERO,YMIN)
      YMAX=MIN(ONE,YMAX)
      ZMIN=MAX(ZERO,ZMIN)
      ZMAX=MIN(ONE,ZMAX)
      CALL SENEXT
C
      IF (.NOT.QCOOR) THEN
      CALL WRNDIE(-5,'SELRPN','no coordinates present')
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
      ELSE
      DO IAT=1,NATOM
C
      QSTACK(QSTART+IAT-1)=.FALSE.
C
C loop over non-crystallographic symmetry operators
      DO INSYM=1,XNNSYM
C
C compute non-crystallographic symmetry related atom (orthogonal coords)
      CALL NCSXYZ(1,XC(IAT),YC(IAT),ZC(IAT),INSYM,XRRF,YRRF,ZRRF)
C
C loop over crystallographic symmetry operators
      DO ISYM=1,XRNSYM
C
C compute symmetry related molecule (in fractional coordinates)
      CALL XPSYMO(1,XRRF,YRRF,ZRRF,ISYM,XRRFF,YRRFF,ZRRFF)
C
C apply mod operation
      XRRFF=MOD(XRRFF+BIG,ONE)
      YRRFF=MOD(YRRFF+BIG,ONE)
      ZRRFF=MOD(ZRRFF+BIG,ONE)
      IF          (XRRFF.GE.XMIN.AND.XRRFF.LE.XMAX
     &        .AND.YRRFF.GE.YMIN.AND.YRRFF.LE.YMAX
     &        .AND.ZRRFF.GE.ZMIN.AND.ZRRFF.LE.ZMAX) THEN
      QSTACK(QSTART+IAT-1)=.TRUE.
      END IF
      END DO
      END DO
      END DO
      END IF
C==================================
      ELSE
      WRITE(6,'(2A)') ' %SELRPN-ERR: Unrecognized word ',WD(1:WDLEN)
      CALL WRNDIE(-5,'SELRPN','parsing error')
C exit selection subroutine
      GOTO 99999
      END IF
C-end-parse-token
C
77777 CONTINUE
C=====================================
      IF (WD(1:4).EQ.'AROU') THEN
C
      IF (.NOT.QCOOR) THEN
      CALL WRNDIE(-5,'SELRPN','no coordinates present')
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
      CALL SENEXT
      ELSE
      CALL NEXTF('AROUnd=',RMAX)
      RMAX2=RMAX*RMAX
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1)=QSTACK(INDX1-QSTKLN)
      END DO
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1-QSTKLN)=.FALSE.
      END DO
C
      QINDX2=QSTART
      QKNOWN=.TRUE.
      DO IAT2=1,NATOM
      IF (QSTACK(QINDX2)) THEN
      IF (.NOT. INITIA(IAT2,X,Y,Z)) THEN
      QKNOWN=.FALSE.
      ELSE
      QINDX1=QSTART-QSTKLN-1
C
C================
C-VECTORIZED-CODE
C================
      DO IAT=1,NATOM
      DXL=XC(IAT)-XC(IAT2)
      DYL=YC(IAT)-YC(IAT2)
      DZL=ZC(IAT)-ZC(IAT2)
      DIST2=DXL*DXL+DYL*DYL+DZL*DZL
      QINDX1=QINDX1+1
      QSTACK(QINDX1)=QSTACK(QINDX1).OR.(DIST2.LE.RMAX2)
      END DO
      END IF
      END IF
      QINDX2=QINDX2+1
      END DO
      END IF
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
      IF (.NOT.QKNOWN) THEN
      WRITE(6,'(A)')
     & ' %SELRPN-AROUnd-ERR: some atom coordinates unknown'
      END IF
      CALL SENEXT
C=====================================
      ELSE IF (WD(1:4).EQ.'SARO') THEN
C
      IF (.NOT.QCOOR) THEN
      CALL WRNDIE(-5,'SELRPN','no coordinates present')
      DO IAT=1,NATOM
      QSTACK(QSTART+IAT-1)=.FALSE.
      END DO
      CALL SENEXT
      ELSE
      CALL NEXTF('SAROund=',RMAX)
      RMAX2=RMAX*RMAX
      CALL SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1)=QSTACK(INDX1-QSTKLN)
      END DO
      DO INDX1=QSTART,QSTART+QSTKLN-1
      QSTACK(INDX1-QSTKLN)=.FALSE.
      END DO
C
C allocate temporary space
      XF=ALLHP(IREAL8(NATOM))
      YF=ALLHP(IREAL8(NATOM))
      ZF=ALLHP(IREAL8(NATOM))
      XFF=ALLHP(IREAL8(NATOM))
      YFF=ALLHP(IREAL8(NATOM))
      ZFF=ALLHP(IREAL8(NATOM))
      XFFF=ALLHP(IREAL8(NATOM))
      YFFF=ALLHP(IREAL8(NATOM))
      ZFFF=ALLHP(IREAL8(NATOM))
C
C loop over non-crystallographic symmetry operators
      DO INSYM=1,XNNSYM
C
C compute non-crystallographic symmetry related molecule (orthogonal coords)
      CALL NCSXYZ(NATOM,X,Y,Z,INSYM,HEAP(XF),HEAP(YF),HEAP(ZF))
C
C loop over crystallographic symmetry operators
      DO ISYM=1,XRNSYM
C
C compute symmetry related molecule (in fractional coordinates)
      CALL XPSYMO(NATOM,HEAP(XF),HEAP(YF),HEAP(ZF),ISYM,
     &            HEAP(XFF),HEAP(YFF),HEAP(ZFF))
C
      QINDX2=QSTART
      QKNOWN=.TRUE.
      DO IAT2=1,NATOM
      IF (QSTACK(QINDX2)) THEN
      IF (.NOT. INITIA(IAT2,X,Y,Z)) THEN
      QKNOWN=.FALSE.
      ELSE
      QINDX1=QSTART-QSTKLN-1
C
C compute the minimum image distance
      CALL XPIMAG(X(IAT2),Y(IAT2),Z(IAT2),
     &            1,NATOM,HEAP(XFF),HEAP(YFF),HEAP(ZFF),
     &            HEAP(XFFF),HEAP(YFFF),HEAP(ZFFF),.FALSE.,0)
      CALL SELMMM(NATOM,HEAP(XFFF),HEAP(YFFF),HEAP(ZFFF),QSTACK,
     &            QINDX1,RMAX2)
      END IF
      END IF
      QINDX2=QINDX2+1
      END DO
C
      END DO
      END DO
C
C deallocate temporary space
      CALL FREHP(XF,IREAL8(NATOM))
      CALL FREHP(YF,IREAL8(NATOM))
      CALL FREHP(ZF,IREAL8(NATOM))
      CALL FREHP(XFF,IREAL8(NATOM))
      CALL FREHP(YFF,IREAL8(NATOM))
      CALL FREHP(ZFF,IREAL8(NATOM))
      CALL FREHP(XFFF,IREAL8(NATOM))
      CALL FREHP(YFFF,IREAL8(NATOM))
      CALL FREHP(ZFFF,IREAL8(NATOM))
C
      END IF
      CALL SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
      IF (.NOT.QKNOWN) THEN
      WRITE(6,'(A)')
     & ' %SELRPN-SAROund-ERR: some atom coordinates unknown'
      END IF
C
      CALL SENEXT
C=====================================
      END IF
C
      CALL SEUP(LEVEL,RETMAX)
C     return to address:
      IF (RETADR(LEVEL).EQ.5) GOTO 3335
      IF (RETADR(LEVEL).EQ.6) GOTO 3336
      IF (RETADR(LEVEL).EQ.7) GOTO 3337
      IF (RETADR(LEVEL).EQ.8) GOTO 3338
      IF (RETADR(LEVEL).EQ.9) GOTO 3339
      IF (RETADR(LEVEL).EQ.10) GOTO 3340
      WRITE(6,'(A)')
     &' %SELRPN-ERR: Return address unknown. Fatal coding error.'
      CALL DIE
C---- END PROCEDURE FACTOR ----------------------------------------------
C
C
C
C
9999  CONTINUE
      QINDX1=QSTART
      DO IAT=1,NATOM
      IF (QSTACK(QINDX1)) THEN
      FLAGS(IAT)=1
      ELSE
      FLAGS(IAT)=0
      END IF
      QINDX1=QINDX1+1
      END DO
C
C error exit label
99999 CONTINUE
      NSELCT=0
      DO IAT=1,NATOM
      IF (FLAGS(IAT).EQ.1) NSELCT=NSELCT+1
      END DO
      IF (WRNLEV.GE.5) THEN
      WRITE(6,'(A,I7,A,I7)')
     & ' SELRPN:',NSELCT,' atoms have been selected out of',NATOM
C
      END IF
      DPVAL = NSELCT
      CALL DECLAR( 'SELECT', 'DP', ' ', DCVAL, DPVAL )
C
      RETURN
      END
C
      SUBROUTINE SEDOWN(LEVEL,RETMAX)
C
C go down one level in SELRPN
C
C Author: Axel T. Brunger
C
      IMPLICIT NONE
C I/O
      INTEGER LEVEL, RETMAX
C begin
      IF (LEVEL.GE.RETMAX) THEN
      WRITE(6,'(A)') ' %SELRPN-ERR: Return address overflow (RETMAX)'
      ELSE
      LEVEL=LEVEL+1
      END IF
      RETURN
      END
C===================================================================
      SUBROUTINE SEUP(LEVEL,RETMAX)
C go up one level in SELRPN
C
C Author: Axel T. Brunger
C
      IMPLICIT NONE
C I/O
      INTEGER LEVEL, RETMAX
C begin
      LEVEL=LEVEL-1
      IF (LEVEL.LT.0) THEN
      WRITE(6,'(A)') ' %SELRPN-ERR: Return address underflow. Check'
      END IF
      RETURN
      END
C===================================================================
      SUBROUTINE SENEXT
C
C get next word  ( lets be able to type HELP at any place !! )
C Author: Axel T. Brunger
C ======================
      IMPLICIT NONE
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
C begin
      CALL NEXTWD('SELRPN>')
C
      DO WHILE (WD(1:4).EQ.'HELP')
C
      CALL CNSHELP('atom-selection')
C
      CALL NEXTWD('SELRPN>')
      END DO
      RETURN
      END
C=====================================================================
      SUBROUTINE SEPUSH(QLIFT,QSTART,QSTKLN,ISTK,STKSIZ)
C
C push stack in SELRPN
C Author: Axel T. Brunger
C ======================
      IMPLICIT NONE
      INTEGER QLIFT, QSTART, QSTKLN, ISTK, STKSIZ
C begin
      QSTART=QSTART+QSTKLN
      IF (QSTART.GT.STKSIZ) THEN
      CALL WRNDIE(-5,'SELRPN>',
     & 'max. depth of selection expression exceeded')
      ELSE
      QLIFT=QLIFT+1
      ISTK=QLIFT
      END IF
      RETURN
      END
C=====================================================================
      SUBROUTINE SEBACK(QLIFT,QSTART,QSTKLN,ISTK)
C
C back stack in SELRPN
C Author: Axel T. Brunger
C ======================
      IMPLICIT NONE
      INTEGER QLIFT, QSTART, QSTKLN, ISTK
C begin
      QLIFT=QLIFT-1
      QSTART=QSTART-QSTKLN
      IF (QSTART.LE.0) THEN
      CALL WRNDIE(-5,'SELRPN>',
     & 'fatal stack underflow, check code')
      END IF
      RETURN
      END
C=====================================================================
      SUBROUTINE SELATT(ARRAY,QSTART,NATOM,QSTACK,QABS)
C
C routine selects atoms based on comparisons with values in ARRAY
C Author: Axel T. Brunger
C ======================
C
      IMPLICIT NONE
C I/O
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'consta.inc'
      DOUBLE PRECISION ARRAY(*)
      INTEGER QSTART, NATOM
      LOGICAL QSTACK(*), QABS
C local
      INTEGER IROPR, IAT
      DOUBLE PRECISION TOL, CUTOFF, DATAPT
      CHARACTER*2 RELAT
C begin
      IROPR=1
      CALL NEXTQL('operator=')
      RELAT=WD
C
C check if we have a composite relat. operator, such as ">=" or "<=" or "=="
      IF (WD(1:1).EQ.'<'.OR.WD(1:1).EQ.'>') THEN
      CALL NEXTQL('operator=')
      IF (WD(1:1).EQ.'=') THEN
      RELAT(2:2)='='
      ELSE
      CALL SAVEWD
      END IF
      END IF
C
      IF (RELAT.EQ.'LT'.OR.RELAT.EQ.'<') THEN
      IROPR=1
      ELSE IF (RELAT.EQ.'LE'.OR.RELAT.EQ.'<=') THEN
      IROPR=2
      ELSE IF (RELAT.EQ.'GT'.OR.RELAT.EQ.'>') THEN
      IROPR=3
      ELSE IF (RELAT.EQ.'GE'.OR.RELAT.EQ.'>=') THEN
      IROPR=4
      ELSE IF (RELAT.EQ.'EQ'.OR.RELAT.EQ.'=') THEN
      IROPR=5
      ELSE IF (RELAT.EQ.'NE'.OR.RELAT.EQ.'#') THEN
      IROPR=6
      ELSE
      CALL DSPERR('SELRPN','unknown relational operator')
      END IF
      CALL NEXTF('comparison=',CUTOFF)
      TOL=MAX(RSMALL,ABS(RSMALL*CUTOFF))
      DO IAT=1,NATOM
      DATAPT=ARRAY(IAT)
      IF (QABS.AND.DATAPT.LT.0) DATAPT=-DATAPT
      IF ((1).EQ.(IROPR)) THEN
      QSTACK(QSTART+IAT-1)= ((DATAPT-CUTOFF).LT.-TOL)
      ELSE IF ((2).EQ.(IROPR)) THEN
      QSTACK(QSTART+IAT-1)= ((DATAPT-CUTOFF).LE.TOL)
      ELSE IF ((3).EQ.(IROPR)) THEN
      QSTACK(QSTART+IAT-1)= ((DATAPT-CUTOFF).GT.TOL)
      ELSE IF ((4).EQ.(IROPR)) THEN
      QSTACK(QSTART+IAT-1)= ((DATAPT-CUTOFF).GE.-TOL)
      ELSE IF ((5).EQ.(IROPR)) THEN
      QSTACK(QSTART+IAT-1)= (ABS(DATAPT-CUTOFF).LT.TOL)
      ELSE IF ((6).EQ.(IROPR)) THEN
      QSTACK(QSTART+IAT-1)= (ABS(DATAPT-CUTOFF).GE.TOL)
      END IF
      END DO
      RETURN
      END
C
      SUBROUTINE SELRCL(QSTACK,STORE,NATOM)
C
C Subroutine fills QSTACK according to the values in STORE
C QSTACK=(STORE.NE.0)
C
C Author: Axel T. Brunger
C =======================
C
      IMPLICIT NONE
C input/output

      INCLUDE 'consta.inc'
      LOGICAL QSTACK(*)
      DOUBLE PRECISION STORE(*)
      INTEGER NATOM
C local
      INTEGER I
C begin
      DO I=1,NATOM
      QSTACK(I)=(ABS(STORE(I)).GT.RSMALL)
      END DO
      RETURN
      END
C====================================================================
      SUBROUTINE SELTVF(PROMPT,VECTOR)
C
C parses a 3-dimensional vector within selection
C
C For syntax see main parsing loop "HELP"
C
C Author: Axel T. Brunger
C =======================
C =======================
C
      IMPLICIT NONE
C input/output
      INCLUDE 'cns.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'funct.inc'
      DOUBLE PRECISION VECTOR(*)
      CHARACTER*(*) PROMPT
C local
      INTEGER I
      LOGICAL OK
      LOGICAL CLOOP
      CHARACTER*1 SBRA, SKET
      PARAMETER (SBRA='(',SKET=')')
C begin
      ERROR=.FALSE.
C
      CALL NEXTWD(PROMPT)
      IF (WD(1:4).EQ.'=   ') CALL NEXTWD(PROMPT)
      IF (WD(1:4).EQ.'?   ') THEN
      WRITE(6,'(1X,2A,3G12.5,A)') PROMPT,'(',(VECTOR(I),I=1,3),')'
      ELSE
      DO I=1,3
      VECTOR(I)=0.0D0
      END DO
C
      IF (WD(1:1).EQ.SBRA) THEN
      I=0
      CALL NEXTWD(PROMPT)
      CLOOP = .TRUE.
      DO WHILE (CLOOP)
C parse the x,y,z form of <vector>
      IF (I.GE.3) THEN
      ERROR=.TRUE.
      WRITE(6,'(A)') ' %SELTVF-ERR: dimension of vector greater 3'
      ELSE
      I=I+1
      VECTOR(I)=DECODF(WD,WDLEN,OK)
      ERROR=.NOT.OK
      IF (ERROR) THEN
      WRITE(6,'(2A)') ' %SELTVF-ERR: numerical value expected : ',
     1               WD(1:WDLEN)
      END IF
      END IF
      CALL NEXTWD(PROMPT)
      IF (WD(1:1).EQ.SKET.OR.ERROR) CLOOP = .FALSE.
      END DO
C
      ELSE
      WRITE(6,'(A)') ' %SELTVF-ERR: vector starts with a "("'
      END IF
      END IF
      RETURN
      END
C=====================================================================
      SUBROUTINE SELMMM(NATOM,XFF,YFF,ZFF,QSTACK,QINDX1,RMAX2)
C
C Routine checks if distances are within specified cutoff RMAX2
C
C Author: Axel T. Brunger
C
      IMPLICIT NONE
C I/O
      INTEGER NATOM
      DOUBLE PRECISION XFF(*), YFF(*), ZFF(*)
      LOGICAL QSTACK(*)
      INTEGER QINDX1
      DOUBLE PRECISION RMAX2
C local
      INTEGER IAT
      DOUBLE PRECISION DIST2
C begin
      DO IAT=1,NATOM
      DIST2=XFF(IAT)**2+YFF(IAT)**2+ZFF(IAT)**2
      QINDX1=QINDX1+1
      QSTACK(QINDX1)=QSTACK(QINDX1).OR.(DIST2.LE.RMAX2)
      END DO
      RETURN
      END
C=====================================================================
      SUBROUTINE XSCAFIL(OP,NATOM,ARRAY,XRERR)
C
C Routine stores scatter factor OP into atomic array ARRAY.
C
C Author: Axel T. Brunger
C
      IMPLICIT NONE
      INCLUDE 'cns.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'vector.inc'
      INCLUDE 'xcrystal.inc'
      INCLUDE 'xtarget.inc'
C I/O
      CHARACTER*(*) OP
      INTEGER NATOM
      DOUBLE PRECISION ARRAY(*)
      LOGICAL XRERR
C local
      INTEGER INDEX
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
C begin
C
      IF (HPATOF.EQ.0.OR.XRNATF.EQ.0) THEN
      CALL DSPERR ( 'PUTSCA',
     & 'SCATTer database undefined.' )
      XRERR = .TRUE.
      DO INDEX=1,NATOM
      ARRAY(INDEX)=ZERO
      END DO
C
      ELSE
      CALL XSCAFIL2(OP,NATOM,ARRAY,XRERR,XRNATF,HEAP(HPATOF),
     &            HEAP(HPINDF),XRSM,XRSA,XRSB,XRSC,XRFP,XRFDP)
      END IF
      RETURN
      END
C===========================================================
      SUBROUTINE XSCAFIL2(OP,NATOM,ARRAY,XRERR,XRNATF,XRATOF,
     &                  XRINDF,XRSM,XRSA,XRSB,XRSC,XRFP,XRFDP)
C
C
C Author: Axel T. Brunger
C
      IMPLICIT NONE
      INCLUDE 'cns.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'vector.inc'
C I/O
      CHARACTER*(*) OP
      INTEGER NATOM
      DOUBLE PRECISION ARRAY(*)
      LOGICAL XRERR
      INTEGER XRNATF, XRATOF(*), XRINDF(*)
      INTEGER XRSM
      DOUBLE PRECISION XRSA(XRSM,4), XRSB(XRSM,4), XRSC(XRSM)
      DOUBLE PRECISION XRFP(XRSM), XRFDP(XRSM)
C local
      INTEGER IAT, IINDEX, INDEX
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D0)
C begin
C
C loop over all atoms
      DO INDEX=1,NATOM
C we need to lookup INDEX
      IINDEX=0
      DO IAT=1,XRNATF
      IF (XRATOF(IAT).EQ.INDEX) IINDEX=IAT
      END DO
C
      IF ( IINDEX.EQ.0 ) THEN
      ARRAY(INDEX)=ZERO
      ELSEIF (OP.EQ.'A1') THEN
      ARRAY(INDEX)=XRSA(XRINDF(IINDEX),1)
      ELSEIF (OP.EQ.'A2') THEN
      ARRAY(INDEX)=XRSA(XRINDF(IINDEX),2)
      ELSEIF (OP.EQ.'A3') THEN
      ARRAY(INDEX)=XRSA(XRINDF(IINDEX),3)
      ELSEIF (OP.EQ.'A4') THEN
      ARRAY(INDEX)=XRSA(XRINDF(IINDEX),4)
      ELSEIF (OP.EQ.'B1') THEN
      ARRAY(INDEX)=XRSB(XRINDF(IINDEX),1)
      ELSEIF (OP.EQ.'B2') THEN
      ARRAY(INDEX)=XRSB(XRINDF(IINDEX),2)
      ELSEIF (OP.EQ.'B3') THEN
      ARRAY(INDEX)=XRSB(XRINDF(IINDEX),3)
      ELSEIF (OP.EQ.'B4') THEN
      ARRAY(INDEX)=XRSB(XRINDF(IINDEX),4)
      ELSEIF (OP.EQ.'C') THEN
      ARRAY(INDEX)=XRSC(XRINDF(IINDEX))
      ELSEIF (OP.EQ.'FP') THEN
      ARRAY(INDEX)=XRFP(XRINDF(IINDEX))
      ELSEIF (OP.EQ.'FDP') THEN
      ARRAY(INDEX)=XRFDP(XRINDF(IINDEX))
      END IF
C
      END DO
C
      RETURN
      END
C


text:
      subroutine compute_phi2(n_iPHI,n_rPHI,Nx,res)
      implicit none
      integer i
      integer Nx
      real*8 n_iPHI(Nx)
      real*8 n_rPHI(Nx)
      real*8 res(Nx)
      real*8 qb
      do i=1, Nx, 1
      qb = n_rPHI(i) ** 2 + n_iPHI(i) ** 2
      res(i)=qb
      end do
      END

ext: f


text:
module fib_mod
    use, intrinsic :: iso_fortran_env, only : I8 => INT64, DP => REAL64
    implicit none

    private
    public :: fib_recursive, fib_skip_ahead, fib_memoization, &
        fib_iterative, fib_closed_form

contains

    recursive function fib_recursive(n) result(fib)
        implicit none
        integer(Kind=I8), value :: n
        integer(kind=I8) :: fib

        if (n == 0_I8) then
            fib = 0_I8
        else if (n == 1_I8) then
            fib = 1_I8
        else
            fib = fib_recursive(n - 1_I8) + fib_recursive(n - 2_I8)
        end if
    end function fib_recursive

    recursive function fib_skip_ahead(n) result(fib)
        implicit none
        integer(kind=I8), value :: n
        integer(kind=I8) :: fib
        if (n == 0_I8) then
            fib = 0_I8
        else if (n == 1_I8 .or. n == 2_I8) then
            fib = 1_I8
        else if (mod(n, 2_I8) == 0) then
            fib = fib_skip_ahead(n/2_I8)*(fib_skip_ahead(n/2_I8 + 1_I8) + &
                                          fib_skip_ahead(n/2_I8 - 1_I8))
        else
            fib = fib_skip_ahead(n/2_I8 + 1_I8)**2 + fib_skip_ahead(n/2_I8)**2
        end if
    end function fib_skip_ahead

    recursive function fib_memoization(n) result(fib)
        implicit none
        integer(kind=I8), value :: n
        integer(kind=I8) :: fib
        integer(kind=I8), parameter :: buffer_size = 1000000
        integer(kind=I8), dimension(0:buffer_size), save :: buffer = -1_I8

        buffer(0) = 0_I8
        buffer(1) = 1_I8
        if (n <= buffer_size) then
            if (buffer(n) == -1_I8) then
                buffer(n) = fib_memoization(n - 1_I8) + fib_memoization(n - 2_I8)
            end if
            fib = buffer(n)
        else
            fib = fib_memoization(n - 1_I8) + fib_memoization(n - 2_I8)
        end if
    end function fib_memoization

    function fib_iterative(n) result(fib)
        implicit none
        integer(kind=I8), value :: n
        integer(kind=I8) :: fib
        integer(kind=I8) :: fib_n_1, i, tmp

        if (n == 0_I8) then
            fib = 0_I8
        else
            fib_n_1 = 0_I8
            fib = 1_I8
            do i = 2_I8, n
                tmp = fib
                fib = fib + fib_n_1
                fib_n_1 = tmp
            end do
        end if
    end function fib_iterative

    function fib_closed_form(n) result(fib)
        implicit none
        integer(kind=I8), value :: n
        integer(kind=I8) :: fib
        real(kind=DP) :: alpha = (1.0_DP + sqrt(5.0_DP))/2.0_DP, &
                         beta  = (1.0_DP - sqrt(5.0_DP))/2.0_DP, &
                         pre_factor = 1.0_DP/sqrt(5.0_DP)

        fib = nint(pre_factor*(alpha**n - beta**n), kind=I8)
    end function fib_closed_form

end module fib_mod

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
MODULE dfpmin_module

    !
    !   Optimize a function using gradient information using the quasi-Newton
    !   method of Broyden, Fletcher, Goldfarb, and Shanno (BFGS)
    !

    !/* external modules    */

    USE shared_constants

    IMPLICIT NONE

    !/* setup   */

    PRIVATE

    PUBLIC:: dfpmin

CONTAINS
!******************************************************************************
!******************************************************************************
SUBROUTINE dfpmin(func, dfunc, p, gtol, maxiter, stpmx, maxfun, success, message, iter)

    !/* external objects        */

    REAL(our_dble), INTENT(INOUT)   :: p(:)

    REAL(our_dble), INTENT(IN)      :: stpmx
    REAL(our_dble), INTENT(IN)      :: gtol

    INTEGER(our_int), INTENT(IN)    :: maxiter
    INTEGER(our_int), INTENT(IN)    :: maxfun
    INTEGER(our_int), INTENT(OUT)   :: iter

    LOGICAL, INTENT(OUT)            :: success

    CHARACTER(150), INTENT(OUT)     :: message

    INTERFACE

        FUNCTION func(p)

            USE shared_constants

            IMPLICIT NONE

            REAL(our_dble), INTENT(IN)  :: p(:)
            REAL(our_dble)              :: func

        END FUNCTION

        FUNCTION dfunc(p)

            USE shared_constants

            IMPLICIT NONE

            REAL(our_dble), INTENT(IN)  :: p(:)
            REAL(our_dble)              :: dfunc(SIZE(p))

        END FUNCTION

    END INTERFACE

    !/* internal objects        */

    INTEGER(our_int)                :: its
    INTEGER(our_int)                :: i

    REAL(our_dble)                  :: hessin(SIZE(p), SIZE(p))
    REAL(our_dble)                  :: pnew(SIZE(p))
    REAL(our_dble)                  :: hdg(SIZE(p))
    REAL(our_dble)                  :: xi(SIZE(p))
    REAL(our_dble)                  :: dg(SIZE(p))
    REAL(our_dble)                  :: g(SIZE(p))
    REAL(our_dble)                  :: stpmax
    REAL(our_dble)                  :: sumdg
    REAL(our_dble)                  :: sumxi
    REAL(our_dble)                  :: fret
    REAL(our_dble)                  :: den
    REAL(our_dble)                  :: fac
    REAL(our_dble)                  :: fad
    REAL(our_dble)                  :: fae
    REAL(our_dble)                  :: fp

    LOGICAL                         :: check

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    fp = func(p)
    g = dfunc(p)

    ! Initialize Hessian
    hessin = zero_dble
    DO i = 1, SIZE(p)
        hessin(i, i) = one_dble
    END DO

    xi = -g

    stpmax = stpmx * MAX(SQRT(DOT_PRODUCT(p, p)), REAL(size(p), our_dble))

    DO its = 1, maxiter

        iter = its

        CALL lnsrch(p, fp, g, xi, pnew, fret, stpmax, check, func, maxfun)

        IF (num_eval == maxfun) THEN

            success = .False.
            message = 'Maximum number of evaluations reached.'

            RETURN

        END IF

        fp = fret
        xi = pnew - p
        p = pnew

        dg = g
        g = dfunc(p)
        den = max(fret, one_dble)

        IF (MAXVAL(ABS(g) * max(ABS(p), one_dble) / den) < gtol) THEN

            success = .True.
            message = 'Gradient less than requested.'

            RETURN

        END IF

        dg =  g - dg
        hdg = matmul(hessin, dg)
        fac = DOT_PRODUCT(dg, xi)
        fae = DOT_PRODUCT(dg, hdg)
        sumdg = DOT_PRODUCT(dg, dg)
        sumxi = DOT_PRODUCT(xi, xi)

        IF (fac ** 2 > eps * sumdg * sumxi) THEN

            fac = one_dble / fac
            fad = one_dble / fae
            dg = fac * xi - fad * hdg

            hessin = hessin + fac * outerprod(xi, xi) - fad * outerprod(hdg,hdg) + fae * outerprod(dg,dg)

        END IF

        xi = -matmul(hessin,g)

    END DO

    success = .False.
    message = 'Exceeded maximum number of iterations.'

END SUBROUTINE
!******************************************************************************
!******************************************************************************
SUBROUTINE lnsrch(xold, fold, g, p, x, f, stpmax, check, func, maxfun)

    !/* external objects        */

    REAL(our_dble), INTENT(OUT)     :: x(:)
    REAL(our_dble), INTENT(OUT)     :: f

    REAL(our_dble), INTENT(IN)      :: xold(:)
    REAL(our_dble), INTENT(IN)      :: stpmax
    REAL(our_dble), INTENT(IN)      :: fold
    REAL(our_dble), INTENT(IN)      :: g(:)

    REAL(our_dble), INTENT(INOUT)   :: p(:)

    LOGICAL, INTENT(OUT)            :: check

    INTEGER(our_int), INTENT(IN)    :: maxfun

    INTERFACE

        FUNCTION func(x)

        USE shared_constants

        IMPLICIT NONE

        REAL(our_dble)              :: func

        REAL(our_dble), INTENT(IN)  :: x(:)

        END FUNCTION

    END INTERFACE

    !/* internal objects        */

    REAL(our_dble), PARAMETER       :: ALF=1.0e-4_our_dble

    INTEGER(our_int)                :: ndum

    REAL(our_dble)                  :: alamin
    REAL(our_dble)                  :: tmplam
    REAL(our_dble)                  :: slope
    REAL(our_dble)                  :: alam2
    REAL(our_dble)                  :: fold2
    REAL(our_dble)                  :: alam
    REAL(our_dble)                  :: disc
    REAL(our_dble)                  :: pabs
    REAL(our_dble)                  :: rhs1
    REAL(our_dble)                  :: rhs2
    REAL(our_dble)                  :: f2
    REAL(our_dble)                  :: a
    REAL(our_dble)                  :: b

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    ndum = size(g)

    check = .FALSE.

    pabs = SQRT(DOT_PRODUCT(p, p))

    IF (pabs > stpmax) p = p * stpmax / pabs

    slope = DOT_PRODUCT(g, p)

    alamin = eps / MAXVAL(ABS(p)/MAX(ABS(xold), one_dble))

    alam = 1.0

    DO
        x(:) = xold(:) + alam * p(:)

        f = func(x)

        IF (num_eval == maxfun) THEN
            RETURN
        END IF

        IF (alam < alamin) THEN

            x(:) = xold(:)
            check = .True.
            RETURN

        ELSEIF (f <= fold + ALF * alam * slope) THEN

            RETURN

        ELSE

            IF (alam == one_dble) THEN

                tmplam = -slope / (two_dble * (f -fold - slope))

            ELSE

                rhs1 = f - fold - alam * slope
                rhs2 = f2 - fold2 - alam2 * slope
                a = (rhs1 / alam ** 2 - rhs2 / alam2 ** 2) / (alam - alam2)
                b = (-alam2 * rhs1 / alam ** 2 + alam * rhs2 / alam2 ** 2) / (alam - alam2)

                IF (a == zero_dble) THEN

                    tmplam = -slope / (two_dble * b)

                ELSE

                    disc = b * b - three_dble * a * slope

                    IF (disc < zero_dble) THEN
                        disc = 0.001
                    END IF

                    tmplam = (-b + SQRT(disc)) / (three_dble * a)

                END IF

                IF (tmplam > half_dble * alam) tmplam = half_dble * alam

            END IF

        END IF

        alam2 = alam
        f2 = f
        fold2 = fold
        alam = MAX(tmplam, 0.1_our_dble * alam)

    END DO

END SUBROUTINE
!******************************************************************************
!******************************************************************************
PURE FUNCTION outerprod(a,b)

    !/* external objects        */

    REAL(our_dble), INTENT(IN)      :: a(:)
    REAL(our_dble), INTENT(IN)      :: b(:)

    REAL(our_dble)                  :: outerprod(SIZE(a), SIZE(b))

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    outerprod = spread(a, dim=2, ncopies=size(b)) * spread(b, dim=1, ncopies=size(a))

END FUNCTION
!******************************************************************************
!******************************************************************************
END MODULE

max_stars_repo_name: boryana-ilieva/respy

max_issues_repo_name: boryana-ilieva/respy

max_forks_repo_name: boryana-ilieva/respy


text:
      double precision function  DRANG ()
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-04-16 DRANG WVS SQRT(abs(TWO*log(U3))) avoids sqrt(-0.0)
c>> 1994-10-20 DRANG Krogh  Changes to use M77CON
c>> 1994-06-24 DRANG CLL Changed common to use RANC[D/S]1 & RANC[D/S]2.
c>> 1992-03-16 DRANG CLL
c>> 1991-11-26 DRANG CLL Reorganized common. Using RANCM[A/D/S].
c>> 1991-11-22 DRANG CLL Added call to RAN0, and DGFLAG in common.
c>> 1991-01-15 DRANG CLL Reordered common contents for efficiency.
C>> 1990-01-23 DRANG CLL Making names in common same in all subprogams.
C>> 1987-06-09 DRANG CLLawson  Initial code.
C
c     Returns one pseudorandom number from the Gausian (Normal)
c     distribution with zero mean and unit standard deviation.
c     Method taken from Algorithm 334, Comm. ACM, July 1968, p. 498.
c     Implemented at JPL in Univac Fortran V in 1969 by Wiley R. Bunton
c     of JPL and Stephen L. Richie of Heliodyne Corp.
c
c     Adapted to Fortran 77 for the MATH 77 library by C. L. Lawson and
c     S. Y. Chiu, JPL, April 1987, 6/9/87.
C     ------------------------------------------------------------------
c--D replaces "?": ?RANG, ?RANUA, RANC?1, RANC?2, ?PTR, ?NUMS, ?GFLAG
C     RANCD1 and RANCD2 are common blocks.
c     Uses intrinsic functions, log and sqrt.
C     Calls DRANUA to obtain an array of uniform random numbers.
c     Calls RAN0 to initialize DPTR and DGFLAG.
C     ------------------------------------------------------------------
c                        Important common variables
c
c     DPTR [integer] and DGFLAG [logical]
c
c          Will be set to DPTR = 1 and DGFLAG = .false.
c          when RAN0 is called from this subr if this
c          is the first call to RAN0.  Also set to these values when
c          RAN1 or RANPUT is called by a user.
c
c          DGFLAG will be set true on return from this subr when the
c          algorithm has values that it can save to reduce the amount
c          of computation needed the next time this function is
c          referenced.  Will be set false in the contrary case.
c
c          DPTR is the index of the last location used in the
c          common array DNUMS().  This index is counted down.
c
c     DNUMS() [floating point]  Buffer of previously computed uniform
c          random numbers.
C     ------------------------------------------------------------------
      integer M
      parameter(M = 97)
      double precision    ONE, TWO
      parameter(ONE = 1.0D0, TWO = 2.0D0)
      double precision    DNUMS(M), R, S, U3, X, XX, Y, YY
      logical FIRST
      common/RANCD2/DNUMS
      integer DPTR
      logical DGFLAG
      common/RANCD1/DPTR, DGFLAG
      save  /RANCD1/, /RANCD2/, FIRST
      save    R, X, Y
      data  FIRST/.true./
c     ------------------------------------------------------------------
      if(FIRST) then
         FIRST = .false.
         call RAN0
      endif
c
      if (.not. DGFLAG .or. DPTR .eq. 1) then
C
C     Use the Von Neuman rejection method for choosing a random point
c     (X,Y) in the unit circle, X**2 + Y**2 .le. 1.0.
c     Then the angle Theta = arctan(Y/X) is random, uniform in
c     (-Pi/2, Pi/2), and Phi = 2*Theta is random, uniform in (-Pi, Pi).
C     Define S = X**2 + Y**2, then
c     sin(Theta) = Y/sqrt(S),    cos(Theta) = X/sqrt(S),
c     sin(Phi) = 2*X*Y/S,    and cos(Phi) = (X**2 - Y**2)/S.
c
   10    continue
c                              Set X = random, uniform in [0., 1.]
      DPTR = DPTR - 1
      if(DPTR .eq. 0) then
         call DRANUA(DNUMS, M)
         DPTR = M
      endif
            X = DNUMS(DPTR)
c                              Set Y = random, uniform in [-1., 1.]
            DPTR = DPTR - 1
            if(DPTR .eq. 0) then
               call DRANUA(DNUMS,M)
               DPTR = M
            endif
            Y = TWO*DNUMS(DPTR) - ONE
c
            XX=X*X
            YY=Y*Y
            S=XX+YY
         if(S .gt. ONE) go to 10
C
C     Choose R randomly from Chi-squared distribution and
c     normalize with S.
C
c                              Set U3 = random, uniform in [0., 1.]
         DPTR = DPTR - 1
         if(DPTR .eq. 0) then
            call DRANUA(DNUMS,M)
            DPTR = M
         endif
         U3 = DNUMS(DPTR)
c        Changed -TWO*log(U3) to abs(TWO*log(U3)) because Lahey LF90
c        2.00 on a pentium produced -0.0 for -TWO*log(1.0), then got a
c        floating point exception on sqrt(-0.0).
         R = sqrt(abs(TWO*(log(U3))))/S
C
C                                Compute result as  R*Sin(PHI)
C
         DRANG = (XX-YY)*R
         DGFLAG = .true.
         return
      endif
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c        Come here when DGFLAG is true and DPTR is not 1.
c
C                                Compute result as  R*Cos(PHI)
      DRANG = TWO*X*Y*R
      DGFLAG=.false.
      return
      end


text: The Giedt Project is a music scene local band named after the Warren and Leta Giedt Hall.  They have written a song named Giedt, with heavy synth and a sample repeating Giedt as the beat.  (File(Giedt.mp3))

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program main

  use mpi

  integer ( kind = 4 ) error
  integer ( kind = 4 ) id
  integer ( kind = 4 ) p
!
!  Initialize MPI.
!
  call MPI_Init ( error )
!
!  Get the number of processes.
!
  call MPI_Comm_size ( MPI_COMM_WORLD, p, error )
!
!  Get the individual process ID.
!
  call MPI_Comm_rank ( MPI_COMM_WORLD, id, error )
!
!  Every MPI process will print this message.
!
  write ( *, '(a,i1,2x,a)' ) 'P', id, '"Hello, world!"'
!
!  Shut down MPI.
!
  call MPI_Finalize ( error )
end program

max_stars_repo_name: bast/cmake-workshop

max_issues_repo_name: bast/cmake-workshop

max_forks_repo_name: bast/cmake-workshop


text:
: d2*   swap 2* swap 2*c ;
: d2/   2/ swap 2/c swap ;
: -rot  swap >r swap r> ;
: 2swap rot >r rot r> ;
: 2over >r >r 2dup r> r> 2swap ;

: sm/rem  \ d n -- rem quot             \ 6.1.2214
   2dup xor >r  over >r  abs >r dabs r> um/mod
   swap r> 0< if  negate  then
   swap r> 0< if  negate  then ;

: fm/mod  \ d n -- rem quot             \ 6.1.1561
   dup >r  2dup xor >r  dup >r  abs >r dabs r> um/mod
   swap r> 0< if  negate  then
   swap r> 0< if  negate  over if  r@ rot -  swap 1-  then then
   r> drop ;

: (umin)  over over- drop carry ;
: umin   (umin) if swap drop exit then  drop ;
: umax   (umin) if drop exit then  swap drop ;
: /string >r swap r@ + swap r> - ;      \ 17.6.1.0245  a u -- a+1 u-1

: roll                                  \ 2.2980 xu..x0 u -- xu-1..x0 xu
    ?dup if
        1+  frp @  ds>mem  mem>
        1-  over cell -                 ( a u' 'xu )
        dup @ >r  !
        mem>ds drop  r>
    then
;

depth 0 assert
2 3 d2* 6 assert 4 assert
-1 5 d2* 11 assert -2 assert
-5 -7 d2/ -4 assert -3 assert

max_stars_repo_name: uho/chad

max_issues_repo_name: uho/chad

max_forks_repo_name: uho/chad


text:
*----------------------------------------------------------------------*
      subroutine leqc_init(xrsnrm,iroute,
     &       me_opt,me_trv,me_mvp,me_rhs,me_dia,me_met,me_scr,
     &       me_special,nspecial,
     &       nincore,lenbuf,
     &       xbuf1,xbuf2,xbuf3,
     &       flist,depend,use_s,
     &       fspc,nspcfrm,
     &       opti_info,opti_stat,
     &       orb_info,op_info,str_info,strmap_info)
*----------------------------------------------------------------------*
*     initialization of LEQ solver
*----------------------------------------------------------------------*
      implicit none

      include 'stdunit.h'
c      include 'def_filinf.h'
      include 'mdef_operator_info.h'
      include 'def_file_array.h'
      include 'def_optimize_info.h'
      include 'def_optimize_status.h'
      include 'ifc_memman.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_strmapinf.h'
      include 'def_orbinf.h'
      include 'opdim.h'
      include 'def_contraction.h'
      include 'def_formula_item.h'
      include 'def_dependency_info.h'

      integer, parameter ::
     &     ntest = 00

      integer, intent(in) ::
     &     iroute, nincore, lenbuf, nspecial, nspcfrm

      type(me_list_array), intent(inout) ::
     &     me_opt(*), me_dia(*),
     &     me_trv(*), me_mvp(*), me_rhs(*), me_scr(*),
     &     me_special(nspecial)
      type(me_list_array), intent(inout) ::
     &     me_met(*)
c      type(file_array), intent(in) ::
c     &     ffopt(*), fftrv(*), ffmvp(*), ffrhs(*), ffdia(*)

      type(formula_item), intent(inout) ::
     &     flist, fspc(nspcfrm)
      type(dependency_info) ::
     &     depend

      type(optimize_info), intent(in) ::
     &     opti_info
      type(optimize_status), intent(inout), target ::
     &     opti_stat

      type(orbinf), intent(in) ::
     &     orb_info
      type(operator_info), intent(inout) ::
     &     op_info
      type(strinf), intent(in) ::
     &     str_info
      type(strmapinf) ::
     &     strmap_info

      real(8), intent(inout) ::
     &     xrsnrm(opti_info%nroot,opti_info%nopt)
      real(8), intent(inout) ::
     &     xbuf1(*), xbuf2(*), xbuf3(*)

      logical, intent(in) ::
     &     use_s(*)

* local
      logical ::
     &     zero_vec(opti_stat%ndim_vsbsp), trafo
      integer ::
     &     idx, jdx, kdx, iroot, irhs,  nred, nadd, nnew, irecscr,
     &     imet, idamp, nopt, nroot, mxsub, lenmat, job, jopt,
     &     ndim_save, ndel, iopt, lenscr, ifree, restart_mode, nselect
      real(8) ::
     &     cond, xdum, xnrm
      real(8), pointer ::
     &     gred(:), vred(:), mred(:),
     &     xmat1(:), xmat2(:), xvec(:), xret(:)
      integer, pointer ::
     &     ndim_rsbsp, ndim_vsbsp, ndim_ssbsp,
     &     iord_rsbsp(:), iord_vsbsp(:), iord_ssbsp(:),
     &     nwfpar(:),
     &     ipiv(:), iconv(:), idxroot(:), idxselect(:)
      type(file_array), pointer ::
     &     ffvsbsp(:), ffscr(:), ffmet(:)
      type(filinf), target ::
     &     fdum

      integer, external ::
     &     ioptc_get_sbsp_rec
      real(8), external ::
     &     dnrm2

      trafo = .false.

      nopt = opti_info%nopt
      nroot = opti_info%nroot
      mxsub = opti_stat%mxdim_sbsp
c      mred => opti_stat%sbspmat(1:)
c      gred => opti_stat%sbspmat(mxsub**2+1:)
c      vred => opti_stat%sbspmat(2*mxsub**2+1:)
      ndim_rsbsp => opti_stat%ndim_rsbsp
      ndim_vsbsp => opti_stat%ndim_vsbsp
      ndim_ssbsp => opti_stat%ndim_ssbsp
      iord_rsbsp => opti_stat%iord_rsbsp
      iord_vsbsp => opti_stat%iord_vsbsp
      iord_ssbsp => opti_stat%iord_ssbsp
      ffscr => opti_stat%ffscr
      ffvsbsp => opti_stat%ffvsbsp
      nwfpar => opti_info%nwfpar
      zero_vec = .false.

      allocate(ffmet(nopt))

c      if (nopt.gt.1)
c     &     call quit(1,'leqc_init','not yet adapted for nopt>1')
      if (nroot.gt.mxsub)
     &     call quit(1,'leqc_init','?? nroot>mxsub ??')

      if (nincore.ge.2) then

        do iopt = 1, nopt

         !select case(opti_info%typ_prc(iopt))
         if (opti_info%typ_prc(iopt).eq.optinf_prc_file.or.
     &       opti_info%typ_prc(iopt).eq.optinf_prc_mixed.or.
     &       opti_info%typ_prc(iopt).eq.optinf_prc_traf) then
         !case(optinf_prc_file,optinf_prc_mixed,optinf_prc_traf)

          trafo = opti_info%typ_prc(iopt).eq.optinf_prc_traf

          ! read diagonal pre-conditioner
          if (ntest.ge.100) then
            write(lulog,*) 'current ME-list: ',
     &            trim(me_dia(iopt)%mel%label)
          end if
          call vec_from_da(me_dia(iopt)%mel%fhand,1,xbuf2,nwfpar(iopt))
          if (ntest.ge.100)
     &        write(lulog,*) 'xbuf2 norm = ',
     &                       dnrm2(nwfpar(iopt),xbuf2,1)
          do iroot = 1, nroot
            if (trafo) then
              ! the present version works only, if we do this for root
              ! number 1; this is, because we modify the norm below
              ! => we have to switch around the loops
              if (iopt.ne.1)
     &             call quit(1,'leqc_init','route with trafo: problem')
              call optc_traf(me_special(2)%mel,1,xrsnrm(iroot,iopt),
     &                    me_rhs(iopt)%mel,iroot,
     &                    fspc(1),'B',me_special,nspecial,
     &                    nwfpar(iopt),xbuf1,
     &                    orb_info,op_info,str_info,strmap_info)
              call vec_from_da(me_special(2)%mel%fhand,1,xbuf1,
     &                       nwfpar(iopt))
            else
              call vec_from_da(me_rhs(iopt)%mel%fhand,iroot,xbuf1,
     &                       nwfpar(iopt))
            end if
            ! divide rhs's by preconditioner
            if (ntest.ge.100)
     &           write(lulog,*) 'xbuf1 norm = ',
     &                          dnrm2(nwfpar(iopt),xbuf1,1)
c            xnrm = dnrm2(nwfpar(iopt),xbuf1,1)
c            xrsnrm(iroot,iopt) = xnrm
            xnrm = 0d0
            do jopt = 1, nopt
              xnrm = xnrm+xrsnrm(iroot,jopt)**2
            end do
            xnrm = sqrt(xnrm)
            ! %shift is for shifted LEQ
C            call diavc(xbuf1,xbuf1,1d0/xnrm,xbuf2,opti_info%shift,
            call diavc(xbuf1,xbuf1,1d0,xbuf2,opti_info%shift,
     &                 nwfpar(iopt))
            if (ntest.ge.100)
     &           write(lulog,*) 'xbuf1 after division: ' //
     &                          ' norm = ', dnrm2(nwfpar(iopt),xbuf1,1)
            if (trafo) then
              call vec_to_da(me_special(2)%mel%fhand,1,
     &                       xbuf1,nwfpar(iopt))
              call optc_traf(me_opt(iopt)%mel,iroot,xdum,
     &                    me_special(2)%mel,1,
     &                    fspc(1),'F',me_special,nspecial,
     &                    nwfpar(iopt),xbuf1,
     &                    orb_info,op_info,str_info,strmap_info)
              ! copy to scr list
              ! original list was used to ensure spin symmetry if needed
              call switch_mel_record(me_scr(iopt)%mel,iroot)
              call list_copy(me_opt(iopt)%mel,me_scr(iopt)%mel,.false.)
            else
              call vec_to_da(ffscr(iopt)%fhand,iroot,xbuf1,nwfpar(iopt))
            end if
          end do ! iroot
         !case(optinf_prc_blocked)
         else if (opti_info%typ_prc(iopt).eq.optinf_prc_blocked) then
          if (nincore.lt.3)
     &         call quit(1,'leqc_init',
     &         'I need at least 3 incore vectors (prc_blocked)')
          do iroot = 1, nroot
            call vec_from_da(me_rhs(iopt)%mel%fhand,iroot,xbuf1,
     &                       nwfpar(iopt))
c            xnrm = dnrm2(nwfpar(iopt),xbuf1,1)
c            xrsnrm(iroot,iopt) = xnrm
            xnrm = 0d0
            do jopt = 1, nopt
              xnrm = xnrm+xrsnrm(iroot,jopt)**2
            end do
            xnrm = sqrt(xnrm)
            call dscal(nwfpar(iopt),1d0/xnrm,xbuf1,1)
            call optc_prc_special2(me_rhs(iopt)%mel,me_special,
     &                                                      nspecial,
     &                         me_opt(iopt)%mel%op%name,
     &                         opti_info%shift,
     &                         nincore,xbuf1,xbuf2,xbuf3,lenbuf,
     &                         orb_info,op_info,str_info,strmap_info)
            call vec_to_da(ffscr(iopt)%fhand,iroot,xbuf1,nwfpar(iopt))
          end do
         !case default
         else
           call quit(1,'leqc_init','unknown preconditioner type')
         end if
         !end select
        end do
      else

        call quit(1,'leqc_init','incore<2: do the programming')

        ! ... something using:
            call da_diavec(ffscr(iopt)%fhand,iroot,1,0d0,
     &                     ffscr(iopt)%fhand,iroot,1,1d0/xnrm,
     &                      me_dia(iopt)%mel%fhand,1,1,0d0,-1d0,
     &                      nwfpar(iopt),xbuf1,xbuf2,lenbuf)

      end if

      if (ndim_vsbsp.ne.0)
     &     call quit(1,'leqc_init','ndim_vsbsp.ne.0 ???')

      do iopt = 1, nopt
        if (use_s(iopt)) then
          ! assign op. with list containing the scratch trial vector
          call assign_me_list(me_scr(iopt)%mel%label,
     &                        me_opt(iopt)%mel%op%name,op_info)

          ! calculate metric * scratch trial vector
          allocate(xret(depend%ntargets),idxselect(depend%ntargets))
          nselect = 0
          call select_formula_target(idxselect,nselect,
     &                me_met(iopt)%mel%label,depend,op_info)
          do iroot = 1, nroot
            call switch_mel_record(me_met(iopt)%mel,iroot)
            call switch_mel_record(me_scr(iopt)%mel,iroot)
            call frm_sched(xret,flist,depend,idxselect,nselect,
     &             .true.,.false.,op_info,str_info,strmap_info,orb_info)
            me_met(iopt)%mel%fhand%last_mod(iroot) = -1
          end do
          deallocate(xret,idxselect)

          ! reassign op. with list containing trial vector
          call assign_me_list(me_trv(iopt)%mel%label,
     &                        me_opt(iopt)%mel%op%name,op_info)
          ffmet(iopt)%fhand => me_met(iopt)%mel%fhand
        else
          ffmet(iopt)%fhand => fdum
        end if
      end do

      ! orthogonalize initial subspace
      call optc_orthvec(nadd,.false.,
     &                  opti_stat%ffssbsp,iord_ssbsp,1d0, !1d0:dummy
     &                  ffvsbsp,
     &                  iord_vsbsp,ndim_vsbsp,mxsub,zero_vec,
     &                  use_s,0,ffmet,ffscr,nroot,nopt,
     &                  nwfpar,nincore,xbuf1,xbuf2,xbuf3,lenbuf)

      ! set nadd
      if (nadd.eq.0)
     &       call quit(0,'leqc_init',
     &       'solver in problems: only linear dependent '//
     &       'new directions?')
      opti_stat%nadd = nadd

      ! |Mv> subspace organisation should be identical to |v> subsp.
      ndim_rsbsp = ndim_vsbsp
      iord_rsbsp = iord_vsbsp
      ! dto. for |Sv> subspace
      ndim_ssbsp = ndim_vsbsp
      iord_ssbsp = iord_vsbsp

      deallocate(ffmet)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      program t
      integer::ierr
      character(5)::a='*****'
      character(5)::b='*****'
 100  format (i0)
      open(unit=77,file='infile',status='old',iostat=ierr)
      read(77,fmt='(a)',iostat=ierr) a
      print 100,ierr
      print '(a)',a
      rewind(unit=77,iostat=ierr)
      print 100,ierr
      read(77,fmt='(a)',iostat=ierr) b
      print 100,ierr
      print '(a)',b
      close(77,iostat=ierr)
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
## Dear wiki vandal: please dont mess with the phone number.

Americana Arms Apartments are now called Academy Lane Apartments.

Just Renovated! Woodlike vinyl floors, Plush Carpet, New Cabinet Fronts, Refinished Countertops and Bathrooms and New insullation to increase energy efficiency and reduce sound!!

Americana Arms features newly upgraded units. They have 157 onebedroom apartments and 1 huge two bedroom apartment. The deposit is $325.00, and the application is $42.00 per person.  Most people seem to like them, and the management is responsive.  The apartments are large for the price, but have fairly thin walls.  This may not be a problem for you if they have the older, more grad student heavy population that theyve had in the past.

Fiber optic wifi wireless internet provided by Instaconnect  (with the problems that usually entails)   and 90+ channels of Comcast Digital Cable Television are included in the rent.  Americana Arms is a smoking smokefree community. Residents are permitted to smoke 20 feet from an exterior door. The apartments lack dishwashers and garbage disposals.

The complex is within walking distance of grocery shopping Davis Food CoOp and the North Davis Nugget.  The Unitrans E Bus Line and Yolobus 42B and 43 stop right in front of the complex.  For close shopping, Woodland Walmart is an approximate 10 mile drive, and the Dixon Walmart Supercenter is just under 10.  For IKEA, head on over to West Sacramento (13 miles).

 Photos

Older Reviews

   /2009 Reviews Reviews from 2005 to 2009

Current Reviews

20100806 09:22:17 nbsp Ive been very happy here!  The management and maintenance has been awesome!  Any problem minor or major has been attended to very quickly.  A couple of perks: The office will take packages for you so you dont have to try and chase the UPS driver.  They also sell stamps and have a fax machine, which can be convenient. I live on the inside row, so I do not hear any street traffic and the train noise (best of my memory) does not happen in the middle of the night.  I think its only a couple times during the day  hasnt been a problem for me. In fact, I kind of like it.

The pool is great! Users/toolej

20100909 21:18:29 nbsp Weve lived at Americana Arms for the past 3 years and have really enjoyed our stay. The management is very friendly and responsible and maintenance has been absolutely great. Weve lived here for so long, because we think it is the best deal in Davis for a 1bedroom. Wireless and cable are included in the rent, and the wireless is awesome. We can stream our Netflix anytime and weve had no problems whatsoever. The rooms are of decent size and are pretty basic, but there is a nicesized closet in the bedroom. The management recently put in doublepaned windows too, which has helped us save money with energy.

The grounds are well kept, the pool is nice and the location is great since campus is less than 10 minutes away by bike and you can walk downtown. Mostly grad students live here as well, so it doesnt get too loud in the evenings or on the weekends. Overall, weve had almost no complaints with anything. Weve had many friends come and stay with us and theyve all been impressed with our standard of living (especially given the price).

We will be the first to admit that there are nicer places in Davis, but the problem is that youre going to easily pay 23 hundred dollars more per month (2k3k more per year) for something that isnt much of an upgrade. This place has everything you need and more, and is especially great for graduate students and the like.  Users/NickandAdelina

20101011 17:33:32 nbsp Im surprised by all the negative comments left by others. Ive lived here for 2 years and have had a pretty good experience. The apartment is kept in good condition and while it is old, maintenance is great and management is nice. You get a pretty decent space for about $800/month. There arent many other one bedroom apartments in Davis that you can rent for that price. Overall I think this is a nice place! Users/lawrenceucd

20101220 18:13:31 nbsp Note: Since I’m posting from home, my review will apparently show up with from the same IP address that management has since we’re all on one network. I’m too lazy right now to go somewhere else that has wifi to post this, but if anyone actually questions the validity of this review coming from a resident, then I will repost it from school or something. Cheers.

My husband and I have lived in Americana Arms for over half a year now, and since my time here has included experiences of how the apartments feel during different seasons, I feel that I can accurately give an honest and fairly knowledgeable review. One of the reasons I want to leave a review is because all the bad reviews on here almost scared me away from living here.

Overall, I think Americana Arms are very good apartments. They’re cheap, they’re pretty centrally located, the staff and maintenance are nice, and the neighbors are usually easy to get along with. The staff here is professional and has a good attitude. I know that there is a summer intern person, but her name escapes me, but there is Alicia and Leslie in the office, and Roy and some others as maintenance (sorry that I don’t know the other people’s names, I just haven’t interacted with them). Alicia’s really nice, upbeat, and friendly. I think she’s the more peopleperson of the two. Leslie is professional and nice as well, but I never really know how to talk to her and find her a little intimidating. It might just be me, but it always struck me as Alicia is the one who talks to new and prospective residents, and Leslie is the stricter one who would deal with any problems that might arise. This is just an opinion, as I’ve never had any problems during my stay here. Roy’s quiet but nice, and just kind of gets the job done is on his way to the next job to do. Something else I really like about the staff is I feel that they treat me like an adult. When I was searching for apartments, I checked out College Square Apartments and they had this whole “So, when are Mom and Dad going to sign the lease?” attitude that really annoyed me. Americana Arms is geared towards graduate students and other responsible people that just want to live in peace, and as a person who is paying my way through school, I appreciate the respect I get from them.

The apartments are located right next to the train tracks, and you can hear the trains passing by quite clearly. I think management said they only pass by a couple times a day or something like that, but I’ve hear them pass by several times. The noise really isn’t a problem though, and after a while you don’t even notice it anymore, but I suppose it could bother some people. I grew up near some train tracks, so I actually like the sound when I’m paying attention to it. A friend of mine who also lives here has an apartment that faces the tracks, and she likes looking as the train cars pass by, and no one I’ve talked to minds the sound.

The cost of utilities hasn’t been a problem here at all. My husband and I moved to Davis so I could attend UC Davis, and the summer here kicked my butt – at first. Surprising after the few days I spent adjusting to the weather of the valley, we rarely used the A/C. I believe that this is because our apartment windows face north, so the sunlight never directly hits us. We also haven’t really used the heater this winter, except in November when I was sick. All in all, PG&E has been very manageable. However, I find that we’ve used the A/C and heater much less than other residents, so obviously this will depend on what kind of temperature your body needs.

Internet and cable are included in the rent. The internet has been pretty reliable; it seemed faster in the summer before all the students moved in, but even now, it’s fine. The internet has gone out on me for short periods of time, like a few minutes usually, but there were two occasions where it went out for a whole day because something had to be fixed in the office. The cable was kind of hard to set up, but once that was done, we haven’t had any real problems. My neighbors complain that their cable goes out all the time, but the only thing we’ve experienced is a few seconds of the screen going all staticy and then returning to normal.

Communication between residents and staff is prompt and effective. There was only one occasion where my email went unresponded when, during the winter, our pilot went out once during the weekend. There is an emergency afterhours number which I could have called, but since I was dealing with the cold just fine with all my blankets and my significant other beside me and wasn’t really using the heater much anyway, I just sent an email. The email got lost in the system or something, so I sent another on Tuesday (I think), and the pilot was either lit later that same day or the next day.  My first email about the pilot was delivered to the staff’s inbox weeks later, and Alicia actually called me when it did, worried that my pilot had gone unlit for so long. I thought that was very nice of her to check in like that.

For their price, these apartments have a decent amount of space. There’s a walkin closet, but even a person with a moderate amount of belongings can fill up an apartment here pretty quickly. It would be nice if there were outside storage of some sort, but these are old buildings, so some things are to be expected.  Since the building is so old, there are some other things that you should expect when moving here. The walls are thin, so you can hear your neighbors pretty clearly – this usually isn’t a problem as most people make the most noise in the day and quiet down at night to sleep. I’ve only had a problem with noise twice, and both times it was more so the fact that my husband or myself were trying to sleep a little earlier than normal due to work or a cold or something rather than unusual noise at odd hours. There is no garbage disposal, as other people have said, but that shouldn’t be a problem as there’s this little drainercatch thing to strain stuff before it gets down the pipe. There is no dish washer, but I don’t see why that bothers people so much, personally.

Really, American Arms has pros and cons like any place you can live. However, they have a good price, nice staff, the space, and the inclusion of cable and internet in the rent. You can probably find a nicer place in Davis with more space, but overall, I think these are very nice apartments and a good choice to make your home.
 Users/Michellaneous

20110406 10:18:46 nbsp I just moved out of Americana Arms to leave Davis, and I was very pleased with my experience there.  I found from the time of my initial inquiry about renting there, through the leasing process, during my time in the complex, and through my move out, the building manager and office staff were wonderful.  Friendly, professional, flexible, and thoughtful peopleevery time I went to the office, they greeted me by name and helped me with anything I needed.

While living there, I appreciated the incredible responsiveness of the maintenance staff (including a few times late at night when I blew a fuse and needed someone to come trip the breaker for meat 2am, someone was there within 20 minutes!)  The building is clean and wellmaintained, rent is very affordable, neighbors are friendly, and the amenities are greatand of course the location cant be beat!  This is a great place to live.   Users/AshleyHoran

20110902 22:11:58 nbsp Do not live here. I live on the first floor and the walls/floors are paper thin. You can hear the upstairs residents as the floor boards crack ever step they make. Also the internet seems to have severe issues. The server seems to disconnect users ofter and when I say ofter its like every 5 minutes. Connection speads are ridiculously slow if you are lucky to get connected. Not worth the $815 a month I pay. Look elsewhere...you have been warned! Users/davistemp
     Under new management now, who work for a huge heartless corporation, I hear. Didnt even give the previous managers the option to continue working under the new regime, after theyd been working there for years: a quick but not nice sweep...wonder how that attitude will translate to new management/tenant relations?Users/RaoulDuke

I would like to speak in regard to this comment, I believe the management is bighearted and the corporation has nothing to do with how we work on a personal level with our residents, and neighbors. During most property sales or changes in managment, it is pretty standard that the staff change with it. However  most of previous staff  were given options and quit by the way. Several of which were reoffered positions after that as well.
 Users/AngeliqueDonovan

20120826 13:28:25 nbsp Under new management, as RaoulDuke said, and so far its a reign of terror. They routinely give less than 24 hours notice before entering my apartment (for example, by leaving a note at the end of their business day saying they will enter the following morning). Twice that Im aware of they have come without giving any notice whatsoever. Now weve received notice that theyll be coming in ANY TIME IN THE NEXT SIX DAYS to do work and that we need to have all of our furniture moved 2 feet from the wall. For six days. All of this has taken place over a period of less than three months (since they acquired the business). Who knows what else is in store? Users/Psittacosaurus
 You might want to mention California Civil Code section 1954 to your landlord.  Not giving 24+ hours notice prior to entry is a big nono. Users/TomGarberson

Id like to speak on behalf of the Management, the notice to enter as required by law is always posted, and all the construction that is going on is to complete NECESSARY repairs that were ordered during the inspections upon the sale of the property. Any notice that was for a span of time was a courtesy notice to allow you to be warned  in advance, prior to the shorter notice and not intended to be the actual 24 or 48 hr notice to enter that we would post asas required by law. If any entry was made without notice it was due to an emergency or a misunderstanding. The work was not started until the end of July, FAR from 3 months by the way  and the time it takes to complete is due to having to have inspections from the county and city in between each phase of the construction.
Im very sorry you have to deal with it, and I completely understand your frustration, however I believe your exagerative rants are a bit extreme. If you have complaints, I would think they would be better dealt with my making them to the management personally.   Users/AngeliqueDonovan

20130205 14:47:28 nbsp Ive been watching Craigslist. Looks like theyll be changing the name from Americana Arms to Academy Lane? Anyone have more info on this. Used to live here. Looks like with all the upgrades that the rent will also be going UP! Users/jenb
  Yes, they are now Academy Lane Apartments.

20140328 16:01:49 nbsp I have lived here almost four years and really enjoy these clean, affordable apartments.  They are centrally located near downtown Davis and close to Community park, the library, etc.  I bike a lot and it was nice to see the bike racks get modernized a year ago.  In the summer, I make frequent use of the pool and BBQs outside.  Management is very attentive to the needs of the tenants; whenever I have needed something fixed, they take care of it with 24 hours. I also really enjoy the activities they have for tenants around the holiday, such as BBQs, ice cream social, and raffles. Users/DavidMcHugh