C=================================================================
C========================== PMB BNCS =============================
C=================================================================
CMAW	Copywrite 2002, Mark A White, UTMB Galveston
C=================================================================
C======================= HOW TO INSTALL ==========================
C=================================================================
CMAW Copy this file to $CNS_SOLVE/source
CMAW Edit the following lines of xropti.f
CMAW      CALL BNCS2(RENER,TEMP,NATOM,WMAIN,IMOVE,FLAG,XRATOF,
CMAW To read (without the CMAW):
CMAW      CALL BNCS4(RENER,TEMP,NATOM,WMAIN,IMOVE,FLAG,XRATOF,
CMAW and ...
CMAW      CALL XROPT2(HEAP(FLAGS),HEAP(IBSIGM),HEAP(IASIGM),
CMAW To read (without the CMAW):
CMAW      CALL XROPT4(HEAP(FLAGS),HEAP(IBSIGM),HEAP(IASIGM),
CMAW
CMAW cd to $CNS_SOLVE and update the executable
CAMW Add the following to $CNS_SOLVE/{machine-dir}/source/Makefile
CMAW add pmb_bncs.o to OBJS = 
CMAW Add the following line to # source file dependencies
CMAW pmb_bncs.o: pmb_bncs.f
CMAW Add the following line to # include file dependencies 
CMAW pmb_bncs.o: cns.inc comand.inc consta.inc coord.inc funct.inc heap.inc mtf.inc xcrystal.inc xropti.inc xtarget.inc
CMAW make g77update		(Linux g77 compiler)
CMAW make update		(Most UNIX systems)
C=================================================================
C=================================================================
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
      SUBROUTINE BNCS3(RENER,TEMP,NATOM,WMAIN,IMOVE,FLAG,XRATOF,
     &          XRNATF,NATNCS,NUMEQV,POINTR,SIGBNC,NNCSX,BBAR)
C
C routine computes NCS B-factor restraints for a group
C
C Author: Axel T. Brunger
C Modified by MAW to give ENCS= sum( B(1,i) -s(1,2)*B(2,i) )*SIGBNC(k)
C where s(1,2)=avg( B(1,i)/B(2,i) )
C THIS IS THE BSCALE VERSION of BNCS
C
      IMPLICIT NONE
C I/O
      DOUBLE PRECISION RENER, TEMP(*)
      INTEGER NATOM
      DOUBLE PRECISION WMAIN(*)
      INTEGER IMOVE(*), FLAG(*), XRATOF(*), XRNATF
      INTEGER NATNCS, NUMEQV, POINTR(*)
      DOUBLE PRECISION SIGBNC
      INTEGER NNCSX
      DOUBLE PRECISION BBAR(*)
C local
      INTEGER I, J, EQPOS, IPT
      DOUBLE PRECISION DELB, DB
CMAW B-factor scale array max 50 chains
      DOUBLE PRECISION BSCALE(50),BGROUP(50)
C parameter
      DOUBLE PRECISION ZERO, TWO
      PARAMETER (ZERO=0.0D0, TWO=2.0D0)
C
C define pointer array
C fill atom flags
      DO I=1,NATOM
      FLAG(I)=0
      END DO
      DO I=1,XRNATF
      IF (IMOVE(XRATOF(I)).EQ.0) FLAG(XRATOF(I))=I
      END DO
C
      DO I=1,NATNCS
      BBAR(I)=ZERO
      END DO
C
      DO J=1,NUMEQV
      EQPOS=NATNCS*(J-1)
      DO I=1,NATNCS
      BBAR(I)=BBAR(I) + WMAIN(POINTR(I+EQPOS))
      END DO
      END DO
C
      DO I=1,NATNCS
      BBAR(I)=BBAR(I)/NUMEQV
      END DO
C
CMAW 150 Calculate the scale factors for each equivalent chain
      DO J=1,NUMEQV
      BSCALE(J)=ZERO
      BGROUP(J)=ZERO
      END DO
CMAW 155
      DO J=1,NUMEQV
      EQPOS=NATNCS*(J-1)
      DO I=1,NATNCS
      BSCALE(J)=BSCALE(J)+ BBAR(I)/WMAIN(POINTR(I+EQPOS))/NATNCS 
      BGROUP(J)=BGROUP(J)+ (WMAIN(POINTR(I+EQPOS))-BBAR(I))/NATNCS 
      END DO
      END DO
CMAW 162
CMAW      DO J=1,NUMEQV
CMAW      BSCALE(J)=BSCALE(J)/NUMEQV
CMAW      END DO
C 166
C Compute energy and derivatives.
      NNCSX=0
      DO J=1,NUMEQV
      EQPOS=NATNCS*(J-1)
      DO I=1,NATNCS
      IPT=FLAG(POINTR(I+EQPOS))
      IF (IPT.NE.0) THEN
      NNCSX=NNCSX+1
CMAW 175      DELB=WMAIN(POINTR(I+EQPOS))-BBAR(I)
CMAW*BSCALE*      DELB=BSCALE(J)*WMAIN(POINTR(I+EQPOS))-BBAR(I)
CMAW*BGROUP*      DELB=WMAIN(POINTR(I+EQPOS))-BBAR(I)-BGROUP(J)
      DELB=BSCALE(J)*WMAIN(POINTR(I+EQPOS))-BBAR(I)
      RENER=RENER + DELB**2/SIGBNC**2
      DB=TWO*DELB/SIGBNC**2
      TEMP(IPT)=TEMP(IPT) + DB
      END IF
      END DO
      END DO
CMAW PRINT INFO on NCS
      WRITE(6,*)' BNCS3: An improved NCS B-factor Model'
	  WRITE(6,*)' Mark A. White, SCSB, UTMB, (c) 2002:' 
	  WRITE(6,*)' BNCS3: uses SCALED BFACTORS - BGROUP is FYI ' 
      DO J=1,NUMEQV
      WRITE(6,6661)' BSCALE, BGROUP for selection(',J,') = ', 
     &               BSCALE(J),BGROUP(J)
6661  FORMAT( A,I2,A,F10.4, F8.1)
      END DO
CMAW
      RETURN
      END
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
C=====================================================================
C=================================================================
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
      SUBROUTINE BNCS4(RENER,TEMP,NATOM,WMAIN,IMOVE,FLAG,XRATOF,
     &          XRNATF,NATNCS,NUMEQV,POINTR,SIGBNC,NNCSX,BBAR)
C
C routine computes NCS B-factor restraints for a group
C
C Author: Axel T. Brunger
C Modified by MAW to give ENCS= sum( B(1,i) -s(1,2)*B(2,i) )*SIGBNC(k)
C where s(1,2)=avg( B(1,i)/B(2,i) )
C THIS IS THE *BGROUP* VERSION of BNCS
C
      IMPLICIT NONE
C I/O
      DOUBLE PRECISION RENER, TEMP(*)
      INTEGER NATOM
      DOUBLE PRECISION WMAIN(*)
      INTEGER IMOVE(*), FLAG(*), XRATOF(*), XRNATF
      INTEGER NATNCS, NUMEQV, POINTR(*)
      DOUBLE PRECISION SIGBNC
      INTEGER NNCSX
      DOUBLE PRECISION BBAR(*)
C local
      INTEGER I, J, EQPOS, IPT
      DOUBLE PRECISION DELB, DB
CMAW B-factor scale array max 50 chains
      DOUBLE PRECISION BSCALE(50),BGROUP(50)
C parameter
      DOUBLE PRECISION ZERO, TWO
      PARAMETER (ZERO=0.0D0, TWO=2.0D0)
C
C define pointer array
C fill atom flags
      DO I=1,NATOM
      FLAG(I)=0
      END DO
      DO I=1,XRNATF
      IF (IMOVE(XRATOF(I)).EQ.0) FLAG(XRATOF(I))=I
      END DO
C
      DO I=1,NATNCS
      BBAR(I)=ZERO
      END DO
C
      DO J=1,NUMEQV
      EQPOS=NATNCS*(J-1)
      DO I=1,NATNCS
      BBAR(I)=BBAR(I) + WMAIN(POINTR(I+EQPOS))
      END DO
      END DO
C
      DO I=1,NATNCS
      BBAR(I)=BBAR(I)/NUMEQV
      END DO
C
CMAW 150 Calculate the scale factors for each equivalent chain
      DO J=1,NUMEQV
      BSCALE(J)=ZERO
      BGROUP(J)=ZERO
      END DO
CMAW 155
      DO J=1,NUMEQV
      EQPOS=NATNCS*(J-1)
      DO I=1,NATNCS
      BSCALE(J)=BSCALE(J)+ BBAR(I)/WMAIN(POINTR(I+EQPOS))/NATNCS 
      BGROUP(J)=BGROUP(J)+ (WMAIN(POINTR(I+EQPOS))-BBAR(I))/NATNCS 
      END DO
      END DO
CMAW 162
CMAW      DO J=1,NUMEQV
CMAW      BSCALE(J)=BSCALE(J)/NUMEQV
CMAW      END DO
C 166
C Compute energy and derivatives.
      NNCSX=0
      DO J=1,NUMEQV
      EQPOS=NATNCS*(J-1)
      DO I=1,NATNCS
      IPT=FLAG(POINTR(I+EQPOS))
      IF (IPT.NE.0) THEN
      NNCSX=NNCSX+1
CMAW 175      DELB=WMAIN(POINTR(I+EQPOS))-BBAR(I)
CMAW*BSCALE*      DELB=BSCALE(J)*WMAIN(POINTR(I+EQPOS))-BBAR(I)
CMAW*BGROUP*      DELB=WMAIN(POINTR(I+EQPOS))-BBAR(I)-BGROUP(J)
      DELB=WMAIN(POINTR(I+EQPOS))-BBAR(I)-BGROUP(J)
      RENER=RENER + DELB**2/SIGBNC**2
      DB=TWO*DELB/SIGBNC**2
      TEMP(IPT)=TEMP(IPT) + DB
      END IF
      END DO
      END DO
CMAW PRINT INFO on NCS
      WRITE(6,*)' BNCS4: An improved NCS B-factor Model'
	  WRITE(6,*)' Mark A. White, SCSB, UTMB, (c) 2002:' 
	  WRITE(6,*)' BNCS4: uses GROUP delta BFACTORS - BSCALE is FYI ' 
      DO J=1,NUMEQV
      WRITE(6,6661)' BSCALE, BGROUP for selection(',J,') = ', 
     &               BSCALE(J),BGROUP(J)
6661  FORMAT( A,I2,A,F10.4, F8.1)
      END DO
CMAW
      RETURN
      END
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
C=====================================================================
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
      SUBROUTINE XROPT4(FLAGS,BSIGMA,ASIGMA,BGROUP,AGROUP)
C
C See XROPTI above
C
C Author: Axel T. Brunger
C =======================
C
CMAW Modified to use variable bsigma asigma restraints per Tickle and Moss
CMAW Co-author Mark A. White
CMAW
      IMPLICIT NONE
C I/O
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'xropti.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'consta.inc'
      INCLUDE 'xcrystal.inc'
      INCLUDE 'xtarget.inc'
      INTEGER FLAGS(*)
      DOUBLE PRECISION BSIGMA(*), ASIGMA(*)
      INTEGER BGROUP(*), AGROUP(*)
C local
      INTEGER CGSTEP, IAT, NFLAGS
      DOUBLE COMPLEX DBCOMP
      DOUBLE PRECISION STEP, TOLER, TEMP, AVG_B, VSIG
      LOGICAL DEBUG
      INTEGER NBGRP, NAGRP
C parameters
      DOUBLE PRECISION ZERO, ONE, TWO, R15, ASIG_SCALE
      PARAMETER (ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0, R15=1.5D0,
     &           ASIG_SCALE=1.5D0 )
C begin
C
C default values
      CGSTEP=0
C Initial step size for conjugate gradient minimizer
      STEP=10.0D0
C Tolerance on gradient
      TOLER=0.000001D0
C B-factor restraint weight
      RWGHT=-ONE
      QRESTR=.TRUE.
C bond restraint sigma
      DO IAT=1,NATOM
      BSIGMA(IAT)=R15
      END DO
C angle restraint sigma
      DO IAT=1,NATOM
      ASIGMA(IAT)=TWO
      END DO
C minimum allowed B-factor
      BFMIN=2.0D0
C maximum allowed B-factor
      BFMAX=100.D0
C
C number of groups for bond restraints
      NBGRP=0
C atomic selections of groups for bond restraints
      DO IAT=1,NATOM
      BGROUP(IAT)=0
      END DO
C number of groups for angle restraints
      NAGRP=0
C atomic selections of groups for angle restraints
      DO IAT=1,NATOM
      AGROUP(IAT)=0
      END DO
C
C debugging flag
      DEBUG=.FALSE.
C
C parsing
      CALL PUSEND('XROPTI>')
      DO WHILE (.NOT.DONE)
      CALL NEXTWD('XROPTI>')
      CALL MISCOM('XROPTI>',USED)
      IF (.NOT.USED) THEN
C
      IF (WD(1:4).EQ.'HELP') THEN
C
      CALL CNSHELP('cns-xray-optimize-bfactor')
C
C=====================================================================
      ELSE IF (WD(1:4).EQ.'NSTE') THEN
      CALL NEXTI('Number of steps= ',CGSTEP)
C======================================================================
      ELSE IF (WD(1:4).EQ.'STEP'.OR.WD(1:4).EQ.'DROP') THEN
      CALL NEXTF('DROP=',STEP)
C======================================================================
      ELSE IF (WD(1:4).EQ.'TOLE') THEN
      CALL NEXTF('TOLErance=',TOLER)
C======================================================================CMAW
C======================================================================CMAW
      ELSE IF (WD(1:4).EQ.'BVSI') THEN
      CALL SELCTA(FLAGS,NFLAGS,X,Y,Z,.TRUE.)
      CALL MAKIND(FLAGS,NATOM,NFLAGS)
CMAW Calculate the average B value
      AVG_B=ZERO
      DO IAT=1,NFLAGS
      AVG_B=AVG_B+WMAIN(FLAGS(IAT))/NFLAGS
      END DO
CMAW
CMAW This reads the bvsigma_scale value from the input file
CMAW Good values of BVSIGma_scale are 3-5, with 4 being optimal in many cases.
CMAW
      CALL NEXTF('BVSIGma=',VSIG)
      TEMP=AVG_B**2/VSIG
CMAW Generate the variable SigmaB values based on Tickle & Moss
      DO IAT=1,NFLAGS
CMAW      BSIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
CMAW      ASIGMA(FLAGS(IAT))=ASIG_SCALE * WMAIN(FLAGS(IAT))**2/TEMP
      BSIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
      ASIGMA(FLAGS(IAT))=ASIG_SCALE * BSIGMA(FLAGS(IAT))
      END DO
CMAW PRINT INFO on Bsig
      WRITE(6,*)'C===================================================='	   
      WRITE(6,*)' XROPT4: An improved B-factor restraint Model - BVSIG'	   
	  WRITE(6,*)' XROPT4: Mark A. White, SCSB, UTMB, (c) 2002:' 
	  WRITE(6,*)' XROPT4: bond  sigmaB = B^2 / ',TEMP
	  WRITE(6,*)' XROPT4: angle sigmaB = (',ASIG_SCALE,') B^2 / ',TEMP
	  WRITE(6,*)' XROPT4: average B =  ',AVG_B
	  WRITE(6,*)' XROPT4: Input <sigB> =  ',VSIG
      WRITE(6,*)'C===================================================='	   
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-CMAW
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-CMAW
      IF (NFLAGS.GT.0) THEN
      NBGRP=NBGRP+1
      DO IAT=1,NFLAGS
      BGROUP(FLAGS(IAT))=NBGRP
      END DO
      ELSE
      WRITE(6,'(A)') 'XROPT-WRN: zero atom is selected!'
      END IF
C======================================================================
      ELSE IF (WD(1:4).EQ.'BOND') THEN
      TEMP=ZERO
      CALL NEXTF('Bond B-factor Sigma=',TEMP)
      DO IAT=1,NATOM
      BSIGMA(IAT)=TEMP
      END DO
C======================================================================
      ELSE IF (WD(1:4).EQ.'BSIG') THEN
      CALL SELCTA(FLAGS,NFLAGS,X,Y,Z,.TRUE.)
      CALL MAKIND(FLAGS,NATOM,NFLAGS)
CMAW This reads the bsigma_main value from the input file
      CALL NEXTF('BSIGma=',TEMP)
CMAW      DO IAT=1,NFLAGS
CMAW      BSIGMA(FLAGS(IAT))=TEMP
CMAW      END DO
CMAW Generate the BSIGMA values based on Tickle & Moss
      DO IAT=1,NFLAGS
      BSIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
      END DO
CMAW PRINT INFO on Bsig
      WRITE(6,*)' XROPT4: An improved B-factor restraint Model - BSIG'	   
	  WRITE(6,*)' Mark A. White, SCSB, UTMB, (c) 2002:' 
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
      IF (NFLAGS.GT.0) THEN
      NBGRP=NBGRP+1
      DO IAT=1,NFLAGS
      BGROUP(FLAGS(IAT))=NBGRP
      END DO
      ELSE
      WRITE(6,'(A)') 'XROPT-WRN: zero atom is selected!'
      END IF
C======================================================================
      ELSE IF (WD(1:4).EQ.'ANGL') THEN
      TEMP=ZERO
      DO IAT=1,NATOM
      ASIGMA(IAT)=TEMP
      END DO
C======================================================================
      ELSE IF (WD(1:4).EQ.'ASIG') THEN
      CALL SELCTA(FLAGS,NFLAGS,X,Y,Z,.TRUE.)
      CALL MAKIND(FLAGS,NATOM,NFLAGS)
CMAW THis reads the ASIGMA value from the input file
      CALL NEXTF('ASIGma=',TEMP)
CMAW      DO IAT=1,NFLAGS
CMAW      ASIGMA(FLAGS(IAT))=TEMP
CMAW      END DO
CMAW Generate the ASIGMA values based on Tickle & Moss
      DO IAT=1,NFLAGS
      ASIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
      END DO
CMAW PRINT INFO on Bsig
      WRITE(6,*)' XROPT4: An improved B-factor restraint Model - ASIG'	   
	  WRITE(6,*)' Mark A. White, SCSB, UTMB, (c) 2002:' 
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
      IF (NFLAGS.GT.0) THEN
      NAGRP=NAGRP+1
      DO IAT=1,NFLAGS
      AGROUP(FLAGS(IAT))=NAGRP
      END DO
      ELSE
      WRITE(6,'(A)') 'XROPT-WRN: zero atom is selected!'
      END IF
C======================================================================
      ELSE IF (WD(1:4).EQ.'BFMI'.OR.WD(1:4).EQ.'BMIN') THEN
      CALL NEXTF('BMIN=',BFMIN)
C======================================================================
      ELSE IF (WD(1:4).EQ.'BFMA'.OR.WD(1:4).EQ.'BMAX') THEN
      CALL NEXTF('BMAX=',BFMAX)
C======================================================================
      ELSE IF (WD(1:4).EQ.'RWEI') THEN
      CALL NEXTF('B-factor Restraint Weight=',RWGHT)
      QRESTR=ABS(RWGHT).GT.RSMALL
C======================================================================
      ELSE IF (WD(1:4).EQ.'DEBU') THEN
      CALL NEXTLO('DEBUg=',DEBUG)
C======================================================================
      ELSE IF (WD(1:4).EQ.'?   ') THEN
      WRITE(6,'(2A)') ' ---------------------------',
     & 'B-factor-refinement-parameters----------------------'
      WRITE(6,'(A,I6,A,F12.5,A,F12.5,A,F12.5,A,F12.5)')
     &' | NSTEp=',CGSTEP,' DROP=',STEP,' TOLErance=',TOLER,
     &' BMIN=',BFMIN, ' BMAX=',BFMAX
      IF (QRESTR) THEN
      WRITE(6,'(A)')
     &' | restrained B-factor refinement'
C
      IF (RWGHT.LT.ZERO) THEN
      WRITE(6,'(2A)') ' | Weights for B-factor restraints will be set',
     & ' first time by comparing gradients'
      ELSE
      WRITE(6,'(A,F12.5)')
     & ' | Weight for B-factor restraints: ',RWGHT
      END IF
      ELSE
      WRITE(6,'(A)') ' | B-factor restraints turned off'
      END IF
      IF (NBGRP.GT.0) WRITE(6,'(A,I5,A)')
     & ' | There are ',NBGRP,' groups of bond-restraints'
      IF (NAGRP.GT.0) WRITE(6,'(A,I5,A)')
     & ' | There are ',NAGRP,' groups of angle-restraints'
      WRITE(6,'(2A)') ' ---------------------------',
     & '----------------------------------------------------'
C=====================================================================
      ELSE
      CALL CHKEND('XROPTI>',DONE)
      END IF
      END IF
      END DO
      DONE=.FALSE.
C
C carry out minimization or just do debugging of derivatives
      FWMAIN=ALLHP(IREAL8(NATOM))
      CALL XBREF(CGSTEP,STEP,TOLER,DEBUG,BGROUP,AGROUP,NBGRP,NAGRP,
     &           HEAP(FWMAIN),HEAP(HPATOF))
      CALL FREHP(FWMAIN,IREAL8(NATOM))
C
C put the restraint weight into a variable
      CALL DECLAR('B_RWEIGHT','DP',' ',DBCOMP,RWGHT)
C
      RETURN
      END
C======================================================================
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
C=MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW==MAW=
      SUBROUTINE XROPT3(FLAGS,BSIGMA,ASIGMA,BGROUP,AGROUP)
C
C See XROPTI above
C
C Author: Axel T. Brunger
C =======================
C
CMAW Modified to use variable bsigma asigma restraints per Tickle and Moss
CMAW Co-author Mark A. White
CMAW
      IMPLICIT NONE
C I/O
      INCLUDE 'cns.inc'
      INCLUDE 'comand.inc'
      INCLUDE 'xropti.inc'
      INCLUDE 'mtf.inc'
      INCLUDE 'coord.inc'
      INCLUDE 'heap.inc'
      INCLUDE 'funct.inc'
      INCLUDE 'consta.inc'
      INCLUDE 'xcrystal.inc'
      INCLUDE 'xtarget.inc'
      INTEGER FLAGS(*)
      DOUBLE PRECISION BSIGMA(*), ASIGMA(*)
      INTEGER BGROUP(*), AGROUP(*)
C local
      INTEGER CGSTEP, IAT, NFLAGS
      DOUBLE COMPLEX DBCOMP
      DOUBLE PRECISION STEP, TOLER, TEMP, ASIG_SCALE
      LOGICAL DEBUG
      INTEGER NBGRP, NAGRP
C parameters
      DOUBLE PRECISION ZERO, ONE, TWO, R15, AVG_B
      PARAMETER (ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0, R15=1.5D0,
     &           ASIG_SCALE=1.5D0 )
C begin
C
C default values
      CGSTEP=0
C Initial step size for conjugate gradient minimizer
      STEP=10.0D0
C Tolerance on gradient
      TOLER=0.000001D0
C B-factor restraint weight
      RWGHT=-ONE
      QRESTR=.TRUE.
C bond restraint sigma
      DO IAT=1,NATOM
      BSIGMA(IAT)=R15
      END DO
C angle restraint sigma
      DO IAT=1,NATOM
      ASIGMA(IAT)=TWO
      END DO
C minimum allowed B-factor
      BFMIN=2.0D0
C maximum allowed B-factor
      BFMAX=100.D0
C
C number of groups for bond restraints
      NBGRP=0
C atomic selections of groups for bond restraints
      DO IAT=1,NATOM
      BGROUP(IAT)=0
      END DO
C number of groups for angle restraints
      NAGRP=0
C atomic selections of groups for angle restraints
      DO IAT=1,NATOM
      AGROUP(IAT)=0
      END DO
C
C debugging flag
      DEBUG=.FALSE.
C
C parsing
      CALL PUSEND('XROPTI>')
      DO WHILE (.NOT.DONE)
      CALL NEXTWD('XROPTI>')
      CALL MISCOM('XROPTI>',USED)
      IF (.NOT.USED) THEN
C
      IF (WD(1:4).EQ.'HELP') THEN
C
      CALL CNSHELP('cns-xray-optimize-bfactor')
C
C=====================================================================
      ELSE IF (WD(1:4).EQ.'NSTE') THEN
      CALL NEXTI('Number of steps= ',CGSTEP)
C======================================================================
      ELSE IF (WD(1:4).EQ.'STEP'.OR.WD(1:4).EQ.'DROP') THEN
      CALL NEXTF('DROP=',STEP)
C======================================================================
      ELSE IF (WD(1:4).EQ.'TOLE') THEN
      CALL NEXTF('TOLErance=',TOLER)
C======================================================================CMAW
      ELSE IF (WD(1:4).EQ.'BVSI') THEN
      CALL SELCTA(FLAGS,NFLAGS,X,Y,Z,.TRUE.)
      CALL MAKIND(FLAGS,NATOM,NFLAGS)
CMAW Calculate the average B value
      AVG_B=ZERO
      DO IAT=1,NFLAGS
      AVG_B=AVG_B+WMAIN(FLAGS(IAT))/NFLAGS
      END DO
CMAW This reads the bvsigma_scale value from the input file
      CALL NEXTF('BVSIGma=',TEMP)
      TEMP=TEMP*AVG_B
CMAW Generate the BVSIGMA values based on Tickle & Moss
      DO IAT=1,NFLAGS
CMAW      BSIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
CMAW      ASIGMA(FLAGS(IAT))=ASIG_SCALE * WMAIN(FLAGS(IAT))**2/TEMP
      BSIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
      ASIGMA(FLAGS(IAT))=ASIG_SCALE * BSIGMA(FLAGS(IAT))
      END DO
CMAW PRINT INFO on Bsig
      WRITE(6,*)'C===================================================='	   
      WRITE(6,*)' XROPT3: An improved B-factor restraint Model - BVSIG'	   
	  WRITE(6,*)' XROPT3: Mark A. White, SCSB, UTMB, (c) 2002:' 
	  WRITE(6,*)' XROPT3: bond  sigmaB = B^2 / ',TEMP
	  WRITE(6,*)' XROPT3: angle sigmaB = (',ASIG_SCALE,') B^2 / ',TEMP
	  WRITE(6,*)' XROPT3: average B =  ',AVG_B
	  WRITE(6,*)' XROPT3: Input Scale =  ',TEMP/AVG_B
      WRITE(6,*)'C===================================================='	   
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
      IF (NFLAGS.GT.0) THEN
      NBGRP=NBGRP+1
      DO IAT=1,NFLAGS
      BGROUP(FLAGS(IAT))=NBGRP
      END DO
      ELSE
      WRITE(6,'(A)') 'XROPT-WRN: zero atom is selected!'
      END IF
C======================================================================
      ELSE IF (WD(1:4).EQ.'BOND') THEN
      TEMP=ZERO
      CALL NEXTF('Bond B-factor Sigma=',TEMP)
      DO IAT=1,NATOM
      BSIGMA(IAT)=TEMP
      END DO
C======================================================================
      ELSE IF (WD(1:4).EQ.'BSIG') THEN
      CALL SELCTA(FLAGS,NFLAGS,X,Y,Z,.TRUE.)
      CALL MAKIND(FLAGS,NATOM,NFLAGS)
CMAW This reads the bsigma_main value from the input file
      CALL NEXTF('BSIGma=',TEMP)
CMAW      DO IAT=1,NFLAGS
CMAW      BSIGMA(FLAGS(IAT))=TEMP
CMAW      END DO
CMAW Generate the BSIGMA values based on Tickle & Moss
      DO IAT=1,NFLAGS
      BSIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
      END DO
CMAW PRINT INFO on Bsig
      WRITE(6,*)' XROPT3: An improved B-factor restraint Model - BSIG'	   
	  WRITE(6,*)' Mark A. White, SCSB, UTMB, (c) 2002:' 
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
      IF (NFLAGS.GT.0) THEN
      NBGRP=NBGRP+1
      DO IAT=1,NFLAGS
      BGROUP(FLAGS(IAT))=NBGRP
      END DO
      ELSE
      WRITE(6,'(A)') 'XROPT-WRN: zero atom is selected!'
      END IF
C======================================================================
      ELSE IF (WD(1:4).EQ.'ANGL') THEN
      TEMP=ZERO
      DO IAT=1,NATOM
      ASIGMA(IAT)=TEMP
      END DO
C======================================================================
      ELSE IF (WD(1:4).EQ.'ASIG') THEN
      CALL SELCTA(FLAGS,NFLAGS,X,Y,Z,.TRUE.)
      CALL MAKIND(FLAGS,NATOM,NFLAGS)
CMAW THis reads the ASIGMA value from the input file
      CALL NEXTF('ASIGma=',TEMP)
CMAW      DO IAT=1,NFLAGS
CMAW      ASIGMA(FLAGS(IAT))=TEMP
CMAW      END DO
CMAW Generate the ASIGMA values based on Tickle & Moss
      DO IAT=1,NFLAGS
      ASIGMA(FLAGS(IAT))=WMAIN(FLAGS(IAT))**2/TEMP
      END DO
CMAW PRINT INFO on Bsig
      WRITE(6,*)' XROPT3: An improved B-factor restraint Model - ASIG'	   
	  WRITE(6,*)' Mark A. White, SCSB, UTMB, (c) 2002:' 
CMAW*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-
      IF (NFLAGS.GT.0) THEN
      NAGRP=NAGRP+1
      DO IAT=1,NFLAGS
      AGROUP(FLAGS(IAT))=NAGRP
      END DO
      ELSE
      WRITE(6,'(A)') 'XROPT-WRN: zero atom is selected!'
      END IF
C======================================================================
      ELSE IF (WD(1:4).EQ.'BFMI'.OR.WD(1:4).EQ.'BMIN') THEN
      CALL NEXTF('BMIN=',BFMIN)
C======================================================================
      ELSE IF (WD(1:4).EQ.'BFMA'.OR.WD(1:4).EQ.'BMAX') THEN
      CALL NEXTF('BMAX=',BFMAX)
C======================================================================
      ELSE IF (WD(1:4).EQ.'RWEI') THEN
      CALL NEXTF('B-factor Restraint Weight=',RWGHT)
      QRESTR=ABS(RWGHT).GT.RSMALL
C======================================================================
      ELSE IF (WD(1:4).EQ.'DEBU') THEN
      CALL NEXTLO('DEBUg=',DEBUG)
C======================================================================
      ELSE IF (WD(1:4).EQ.'?   ') THEN
      WRITE(6,'(2A)') ' ---------------------------',
     & 'B-factor-refinement-parameters----------------------'
      WRITE(6,'(A,I6,A,F12.5,A,F12.5,A,F12.5,A,F12.5)')
     &' | NSTEp=',CGSTEP,' DROP=',STEP,' TOLErance=',TOLER,
     &' BMIN=',BFMIN, ' BMAX=',BFMAX
      IF (QRESTR) THEN
      WRITE(6,'(A)')
     &' | restrained B-factor refinement'
C
      IF (RWGHT.LT.ZERO) THEN
      WRITE(6,'(2A)') ' | Weights for B-factor restraints will be set',
     & ' first time by comparing gradients'
      ELSE
      WRITE(6,'(A,F12.5)')
     & ' | Weight for B-factor restraints: ',RWGHT
      END IF
      ELSE
      WRITE(6,'(A)') ' | B-factor restraints turned off'
      END IF
      IF (NBGRP.GT.0) WRITE(6,'(A,I5,A)')
     & ' | There are ',NBGRP,' groups of bond-restraints'
      IF (NAGRP.GT.0) WRITE(6,'(A,I5,A)')
     & ' | There are ',NAGRP,' groups of angle-restraints'
      WRITE(6,'(2A)') ' ---------------------------',
     & '----------------------------------------------------'
C=====================================================================
      ELSE
      CALL CHKEND('XROPTI>',DONE)
      END IF
      END IF
      END DO
      DONE=.FALSE.
C
C carry out minimization or just do debugging of derivatives
      FWMAIN=ALLHP(IREAL8(NATOM))
      CALL XBREF(CGSTEP,STEP,TOLER,DEBUG,BGROUP,AGROUP,NBGRP,NAGRP,
     &           HEAP(FWMAIN),HEAP(HPATOF))
      CALL FREHP(FWMAIN,IREAL8(NATOM))
C
C put the restraint weight into a variable
      CALL DECLAR('B_RWEIGHT','DP',' ',DBCOMP,RWGHT)
C
      RETURN
      END
C======================================================================