source/unres/src_CSA/readrtns_csa.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19 c      if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21 c      if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23 c      if (modecalc.eq.14) then
  24 c         call read_MDpar
  25 c         call read_REMDpar
  26 c      endif
  27 C Read MUCA control parameters if required
  28 c      if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39
  40 C Print restraint information
  41 #ifdef MPI
  42       if (.not. out1file .or. me.eq.king) then
  43 #endif
  44       if (nhpb.gt.nss)
  45      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  46      & " distance constraints have been imposed"
  47       do i=nss+1,nhpb
  48         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  49       enddo
  50 #ifdef MPI
  51       endif
  52 #endif
  53 c      print *,"Processor",myrank," leaves READRTNS"
  54       return
  55       end
  56 C-------------------------------------------------------------------------------
  57       subroutine read_control
  58 C
  59 C Read contorl data
  60 C
  61       implicit real*8 (a-h,o-z)
  62       include 'DIMENSIONS'
  63 #ifdef MP
  64       include 'mpif.h'
  65       logical OKRandom, prng_restart
  66       real*8  r1
  67 #endif
  68       include 'COMMON.IOUNITS'
  69       include 'COMMON.TIME1'
  70 c      include 'COMMON.THREAD'
  71       include 'COMMON.SBRIDGE'
  72       include 'COMMON.CONTROL'
  73       include 'COMMON.MCM'
  74 c      include 'COMMON.MAP'
  75       include 'COMMON.HEADER'
  76 c      include 'COMMON.CSA'
  77       include 'COMMON.CHAIN'
  78 c      include 'COMMON.MUCA'
  79 c      include 'COMMON.MD'
  80       include 'COMMON.FFIELD'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 100       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 101       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 102       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 103       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 104       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 105       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 106       call reada(controlcard,'DRMS',drms,0.1D0)
 107       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 108        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 109        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 110        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 111        write (iout,'(a,f10.1)')'DRMS    = ',drms
 112        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 113        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 114       endif
 115       call readi(controlcard,'NZ_START',nz_start,0)
 116       call readi(controlcard,'NZ_END',nz_end,0)
 117       call readi(controlcard,'IZ_SC',iz_sc,0)
 118       timlim=60.0D0*timlim
 119       safety = 60.0d0*safety
 120       timem=timlim
 121       modecalc=0
 122       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 123       minim=(index(controlcard,'MINIMIZE').gt.0)
 124       dccart=(index(controlcard,'CART').gt.0)
 125       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 126       overlapsc=.not.overlapsc
 127       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 128       searchsc=.not.searchsc
 129       sideadd=(index(controlcard,'SIDEADD').gt.0)
 130       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 131       outpdb=(index(controlcard,'PDBOUT').gt.0)
 132       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 133       pdbref=(index(controlcard,'PDBREF').gt.0)
 134       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 135       indpdb=index(controlcard,'PDBSTART')
 136       extconf=(index(controlcard,'EXTCONF').gt.0)
 137       call readi(controlcard,'IPRINT',iprint,0)
 138       call readi(controlcard,'MAXGEN',maxgen,10000)
 139       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 140       call readi(controlcard,"KDIAG",kdiag,0)
 141       call readi(controlcard,"RESCALE_MODE",rescale_mode,0)
 142       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 143      & write (iout,*) "RESCALE_MODE",rescale_mode
 144       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 145       if (index(controlcard,'REGULAR').gt.0.0D0) then
 146         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 147         modecalc=1
 148         refstr=.true.
 149       endif
 150       if (index(controlcard,'CHECKGRAD').gt.0) then
 151         modecalc=5
 152         if (index(controlcard,'CART').gt.0) then
 153           icheckgrad=1
 154         elseif (index(controlcard,'CARINT').gt.0) then
 155           icheckgrad=2
 156         else
 157           icheckgrad=3
 158         endif
 159       elseif (index(controlcard,'THREAD').gt.0) then
 160         modecalc=2
 161         call readi(controlcard,'THREAD',nthread,0)
 162         if (nthread.gt.0) then
 163           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 164         else
 165           if (fg_rank.eq.0)
 166      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 167           stop 'Error termination in Read_Control.'
 168         endif
 169       else if (index(controlcard,'MCMA').gt.0) then
 170         modecalc=3
 171       else if (index(controlcard,'MCEE').gt.0) then
 172         modecalc=6
 173       else if (index(controlcard,'MULTCONF').gt.0) then
 174         modecalc=4
 175       else if (index(controlcard,'MAP').gt.0) then
 176         modecalc=7
 177         call readi(controlcard,'MAP',nmap,0)
 178       else if (index(controlcard,'CSA').gt.0) then
 179         modecalc=8
 180 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 181 crc
 182 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 183 crc
 184 crc        modecalc=9
 185 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 186 cfmc        modecalc=10
 187       else if (index(controlcard,'SOFTREG').gt.0) then
 188         modecalc=11
 189       else if (index(controlcard,'CHECK_BOND').gt.0) then
 190         modecalc=-1
 191       else if (index(controlcard,'TEST').gt.0) then
 192         modecalc=-2
 193       else if (index(controlcard,'MD').gt.0) then
 194         modecalc=12
 195       else if (index(controlcard,'RE ').gt.0) then
 196         modecalc=14
 197       endif
 198
 199       lmuca=index(controlcard,'MUCA').gt.0
 200       call readi(controlcard,'MUCADYN',mucadyn,0)
 201       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 202       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 203      & then
 204        write (iout,*) 'MUCADYN=',mucadyn
 205        write (iout,*) 'MUCASMOOTH=',muca_smooth
 206       endif
 207
 208       iscode=index(controlcard,'ONE_LETTER')
 209       indphi=index(controlcard,'PHI')
 210       indback=index(controlcard,'BACK')
 211       iranconf=index(controlcard,'RAND_CONF')
 212       i2ndstr=index(controlcard,'USE_SEC_PRED')
 213       gradout=index(controlcard,'GRADOUT').gt.0
 214       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 215
 216       if(me.eq.king.or..not.out1file)
 217      & write (iout,'(2a)') diagmeth(kdiag),
 218      &  ' routine used to diagonalize matrices.'
 219       return
 220       end
 221 c------------------------------------------------------------------------------
 222       subroutine molread
 223 C
 224 C Read molecular data.
 225 C
 226       implicit real*8 (a-h,o-z)
 227       include 'DIMENSIONS'
 228 #ifdef MPI
 229       include 'mpif.h'
 230       integer error_msg
 231 #endif
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.GEO'
 234       include 'COMMON.VAR'
 235       include 'COMMON.INTERACT'
 236       include 'COMMON.LOCAL'
 237       include 'COMMON.NAMES'
 238       include 'COMMON.CHAIN'
 239       include 'COMMON.FFIELD'
 240       include 'COMMON.SBRIDGE'
 241       include 'COMMON.HEADER'
 242       include 'COMMON.CONTROL'
 243 c      include 'COMMON.DBASE'
 244 c      include 'COMMON.THREAD'
 245       include 'COMMON.CONTACTS'
 246       include 'COMMON.TORCNSTR'
 247       include 'COMMON.TIME1'
 248       include 'COMMON.BOUNDS'
 249 c      include 'COMMON.MD'
 250 c      include 'COMMON.REMD'
 251       include 'COMMON.SETUP'
 252       character*4 sequence(maxres)
 253       integer rescode
 254       double precision x(maxvar)
 255       character*256 pdbfile
 256       character*320 weightcard
 257       character*80 weightcard_t,ucase
 258       dimension itype_pdb(maxres)
 259       common /pizda/ itype_pdb
 260       logical seq_comp,fail
 261       double precision energia(0:n_ene)
 262       integer ilen
 263       external ilen
 264 C
 265 C Body
 266 C
 267 C Read weights of the subsequent energy terms.
 268
 269        call card_concat(weightcard)
 270        call reada(weightcard,'WLONG',wlong,1.0D0)
 271        call reada(weightcard,'WSC',wsc,wlong)
 272        call reada(weightcard,'WSCP',wscp,wlong)
 273        call reada(weightcard,'WELEC',welec,1.0D0)
 274        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 275        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 276        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 277        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 278        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 279        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 280        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 281        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 282        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 283        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 284        call reada(weightcard,'WBOND',wbond,1.0D0)
 285        call reada(weightcard,'WTOR',wtor,1.0D0)
 286        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 287        call reada(weightcard,'WANG',wang,1.0D0)
 288        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 289 C     Juyong
 290        call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
 291        call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
 292        call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
 293        call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
 294 C
 295        call reada(weightcard,'SCAL14',scal14,0.4D0)
 296        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 297        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 298        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 299        call reada(weightcard,'TEMP0',temp0,300.0d0)
 300        if (index(weightcard,'SOFT').gt.0) ipot=6
 301 C 12/1/95 Added weight for the multi-body term WCORR
 302        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 303        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 304        weights(1)=wsc
 305        weights(2)=wscp
 306        weights(3)=welec
 307        weights(4)=wcorr
 308        weights(5)=wcorr5
 309        weights(6)=wcorr6
 310        weights(7)=wel_loc
 311        weights(8)=wturn3
 312        weights(9)=wturn4
 313        weights(10)=wturn6
 314        weights(11)=wang
 315        weights(12)=wscloc
 316        weights(13)=wtor
 317        weights(14)=wtor_d
 318        weights(15)=wstrain
 319        weights(16)=wvdwpp
 320        weights(17)=wbond
 321        weights(18)=scal14
 322        weights(21)=wsccor
 323 C     JUYONG
 324        weights(24)=wdfa_dist
 325        weights(25)=wdfa_tor
 326        weights(26)=wdfa_nei
 327        weights(27)=wdfa_beta
 328 C
 329
 330       if(me.eq.king.or..not.out1file)
 331      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 332      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 333      &  wturn4,wturn6,
 334      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 335
 336    10 format (/'Energy-term weights (unscaled):'//
 337      & 'WSCC=   ',f10.6,' (SC-SC)'/
 338      & 'WSCP=   ',f10.6,' (SC-p)'/
 339      & 'WELEC=  ',f10.6,' (p-p electr)'/
 340      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 341      & 'WBOND=  ',f10.6,' (stretching)'/
 342      & 'WANG=   ',f10.6,' (bending)'/
 343      & 'WSCLOC= ',f10.6,' (SC local)'/
 344      & 'WTOR=   ',f10.6,' (torsional)'/
 345      & 'WTORD=  ',f10.6,' (double torsional)'/
 346      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 347      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 348      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 349      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 350      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 351      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 352      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 353      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 354      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 355      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 356      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 357      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 358      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 359
 360       if(me.eq.king.or..not.out1file)then
 361        if (wcorr4.gt.0.0d0) then
 362         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 363      &   'between contact pairs of peptide groups'
 364         write (iout,'(2(a,f5.3/))')
 365      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 366      &  'Range of quenching the correlation terms:',2*delt_corr
 367        else if (wcorr.gt.0.0d0) then
 368         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 369      &   'between contact pairs of peptide groups'
 370        endif
 371        write (iout,'(a,f8.3)')
 372      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 373        write (iout,'(a,f8.3)')
 374      &  'General scaling factor of SC-p interactions:',scalscp
 375       endif
 376       r0_corr=cutoff_corr-delt_corr
 377       do i=1,20
 378         aad(i,1)=scalscp*aad(i,1)
 379         aad(i,2)=scalscp*aad(i,2)
 380         bad(i,1)=scalscp*bad(i,1)
 381         bad(i,2)=scalscp*bad(i,2)
 382       enddo
 383 c      call rescale_weights(t_bath)
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 387      &  wturn4,wturn6,
 388      &  wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
 389
 390    22 format (/'Energy-term weights (scaled):'//
 391      & 'WSCC=   ',f10.6,' (SC-SC)'/
 392      & 'WSCP=   ',f10.6,' (SC-p)'/
 393      & 'WELEC=  ',f10.6,' (p-p electr)'/
 394      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 395      & 'WBOND=  ',f10.6,' (stretching)'/
 396      & 'WANG=   ',f10.6,' (bending)'/
 397      & 'WSCLOC= ',f10.6,' (SC local)'/
 398      & 'WTOR=   ',f10.6,' (torsional)'/
 399      & 'WTORD=  ',f10.6,' (double torsional)'/
 400      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 401      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 402      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 403      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 404      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 405      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 406      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 407      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 408      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 409      & 'WDFA_D= ',f10.6,' (DFA, distance)'   /
 410      & 'WDFA_T= ',f10.6,' (DFA, torsional)'   /
 411      & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)'   /
 412      & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
 413
 414       if(me.eq.king.or..not.out1file)
 415      & write (iout,*) "Reference temperature for weights calculation:",
 416      &  temp0
 417       call reada(weightcard,"D0CM",d0cm,3.78d0)
 418       call reada(weightcard,"AKCM",akcm,15.1d0)
 419       call reada(weightcard,"AKTH",akth,11.0d0)
 420       call reada(weightcard,"AKCT",akct,12.0d0)
 421       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 422       call reada(weightcard,"V2SS",v2ss,7.61d0)
 423       call reada(weightcard,"V3SS",v3ss,13.7d0)
 424       call reada(weightcard,"EBR",ebr,-5.50D0)
 425       if(me.eq.king.or..not.out1file) then
 426        write (iout,*) "Parameters of the SS-bond potential:"
 427        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 428      & " AKCT",akct
 429        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 430        write (iout,*) "EBR",ebr
 431        print *,'indpdb=',indpdb,' pdbref=',pdbref
 432       endif
 433       if (indpdb.gt.0 .or. pdbref) then
 434         read(inp,'(a)') pdbfile
 435         if(me.eq.king.or..not.out1file)
 436      &   write (iout,'(2a)') 'PDB data will be read from file ',
 437      &   pdbfile(:ilen(pdbfile))
 438         open(ipdbin,file=pdbfile,status='old',err=33)
 439         goto 34
 440   33    write (iout,'(a)') 'Error opening PDB file.'
 441         stop
 442   34    continue
 443 c        print *,'Begin reading pdb data'
 444         call readpdb
 445 c        print *,'Finished reading pdb data'
 446         if(me.eq.king.or..not.out1file)
 447      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 448      &   ' nstart_sup=',nstart_sup
 449         do i=1,nres
 450           itype_pdb(i)=itype(i)
 451         enddo
 452         close (ipdbin)
 453         nnt=nstart_sup
 454         nct=nstart_sup+nsup-1
 455         call contact(.false.,ncont_ref,icont_ref,co)
 456
 457         if (sideadd) then
 458          if(me.eq.king.or..not.out1file)
 459      &    write(iout,*)'Adding sidechains'
 460          maxsi=1000
 461          do i=2,nres-1
 462           iti=itype(i)
 463           if (iti.ne.10) then
 464             nsi=0
 465             fail=.true.
 466             do while (fail.and.nsi.le.maxsi)
 467 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 468               nsi=nsi+1
 469             enddo
 470             if(fail) write(iout,*)'Adding sidechain failed for res ',
 471      &              i,' after ',nsi,' trials'
 472           endif
 473          enddo
 474         endif
 475       endif
 476
 477       if (indpdb.eq.0) then
 478 C Read sequence if not taken from the pdb file.
 479         read (inp,*) nres
 480 c        print *,'nres=',nres
 481         if (iscode.gt.0) then
 482           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 483         else
 484           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 485         endif
 486 C Convert sequence to numeric code
 487         do i=1,nres
 488           itype(i)=rescode(i,sequence(i),iscode)
 489         enddo
 490 C Assign initial virtual bond lengths
 491         do i=2,nres
 492           vbld(i)=vbl
 493           vbld_inv(i)=vblinv
 494         enddo
 495         do i=2,nres-1
 496           vbld(i+nres)=dsc(itype(i))
 497           vbld_inv(i+nres)=dsc_inv(itype(i))
 498 c          write (iout,*) "i",i," itype",itype(i),
 499 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 500         enddo
 501       endif
 502 c      print *,nres
 503 c      print '(20i4)',(itype(i),i=1,nres)
 504       do i=1,nres
 505 #ifdef PROCOR
 506         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 507 #else
 508         if (itype(i).eq.21) then
 509 #endif
 510           itel(i)=0
 511 #ifdef PROCOR
 512         else if (itype(i+1).ne.20) then
 513 #else
 514         else if (itype(i).ne.20) then
 515 #endif
 516           itel(i)=1
 517         else
 518           itel(i)=2
 519         endif
 520       enddo
 521       if(me.eq.king.or..not.out1file)then
 522        write (iout,*) "ITEL"
 523        do i=1,nres-1
 524          write (iout,*) i,itype(i),itel(i)
 525        enddo
 526        print *,'Call Read_Bridge.'
 527       endif
 528       call read_bridge
 529 C 8/13/98 Set limits to generating the dihedral angles
 530       do i=1,nres
 531         phibound(1,i)=-pi
 532         phibound(2,i)=pi
 533       enddo
 534       read (inp,*) ndih_constr
 535       if (ndih_constr.gt.0) then
 536         read (inp,*) ftors
 537         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 538         if(me.eq.king.or..not.out1file)then
 539          write (iout,*)
 540      &   'There are',ndih_constr,' constraints on phi angles.'
 541          do i=1,ndih_constr
 542           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 543          enddo
 544         endif
 545         do i=1,ndih_constr
 546           phi0(i)=deg2rad*phi0(i)
 547           drange(i)=deg2rad*drange(i)
 548         enddo
 549         if(me.eq.king.or..not.out1file)
 550      &   write (iout,*) 'FTORS',ftors
 551         do i=1,ndih_constr
 552           ii = idih_constr(i)
 553           phibound(1,ii) = phi0(i)-drange(i)
 554           phibound(2,ii) = phi0(i)+drange(i)
 555         enddo
 556       endif
 557       nnt=1
 558 #ifdef MPI
 559       if (me.eq.king) then
 560 #endif
 561        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 562        do i=1,nres
 563          write (iout,'(a3,i5,2f10.1)')
 564      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 565        enddo
 566 #ifdef MP
 567       endif
 568 #endif
 569       nct=nres
 570 cd      print *,'NNT=',NNT,' NCT=',NCT
 571       if (itype(1).eq.21) nnt=2
 572       if (itype(nres).eq.21) nct=nct-1
 573
 574 C     Juyong:READ init_vars
 575 C     Initialize variables!
 576 C     Juyong:READ read_info
 577 C     READ fragment information!!
 578 C     both routines should be in dfa.F file!!
 579
 580 <<<<<<< HEAD
 581       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 582      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 583 =======
 584       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0
 585      &      .and. wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 586 >>>>>>> adam
 587        call init_dfa_vars
 588        print*, 'init_dfa_vars finished!'
 589        call read_dfa_info
 590        print*, 'read_dfa_info finished!'
 591       endif
 592 C
 593 C
 594
 595
 596       if (pdbref) then
 597         if(me.eq.king.or..not.out1file)
 598      &   write (iout,'(a,i3)') 'nsup=',nsup
 599         nstart_seq=nnt
 600         if (nsup.le.(nct-nnt+1)) then
 601           do i=0,nct-nnt+1-nsup
 602             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 603               nstart_seq=nnt+i
 604               goto 111
 605             endif
 606           enddo
 607           write (iout,'(a)')
 608      &            'Error - sequences to be superposed do not match.'
 609           stop
 610         else
 611           do i=0,nsup-(nct-nnt+1)
 612             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 613      &      then
 614               nstart_sup=nstart_sup+i
 615               nsup=nct-nnt+1
 616               goto 111
 617             endif
 618           enddo
 619           write (iout,'(a)')
 620      &            'Error - sequences to be superposed do not match.'
 621         endif
 622   111   continue
 623         if (nsup.eq.0) nsup=nct-nnt
 624         if (nstart_sup.eq.0) nstart_sup=nnt
 625         if (nstart_seq.eq.0) nstart_seq=nnt
 626         if(me.eq.king.or..not.out1file)
 627      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 628      &                 ' nstart_seq=',nstart_seq
 629       endif
 630 c--- Zscore rms -------
 631       if (nz_start.eq.0) nz_start=nnt
 632       if (nz_end.eq.0 .and. nsup.gt.0) then
 633         nz_end=nnt+nsup-1
 634       else if (nz_end.eq.0) then
 635         nz_end=nct
 636       endif
 637       if(me.eq.king.or..not.out1file)then
 638        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 639        write (iout,*) 'IZ_SC=',iz_sc
 640       endif
 641 c----------------------
 642       call init_int_table
 643       if (refstr) then
 644         if (.not.pdbref) then
 645           call read_angles(inp,*38)
 646           goto 39
 647    38     write (iout,'(a)') 'Error reading reference structure.'
 648 #ifdef MPI
 649           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 650           stop 'Error reading reference structure'
 651 #endif
 652    39     call chainbuild
 653           call setup_var
 654 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 655           nstart_sup=nnt
 656           nstart_seq=nnt
 657           nsup=nct-nnt+1
 658           do i=1,2*nres
 659             do j=1,3
 660               cref(j,i)=c(j,i)
 661             enddo
 662           enddo
 663           call contact(.true.,ncont_ref,icont_ref,co)
 664         endif
 665 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 666         call flush(iout)
 667         if (constr_dist.gt.0) call read_dist_constr
 668 c        write (iout,*) "After read_dist_constr nhpb",nhpb
 669         call hpb_partition
 670         if(me.eq.king.or..not.out1file)
 671      &   write (iout,*) 'Contact order:',co
 672         if (pdbref) then
 673         if(me.eq.king.or..not.out1file)
 674      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 675         do i=1,ncont_ref
 676           do j=1,2
 677             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 678           enddo
 679           if(me.eq.king.or..not.out1file)
 680      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 681      &     icont_ref(1,i),' ',
 682      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 683         enddo
 684         endif
 685       endif
 686       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 687      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 688      &    modecalc.ne.10) then
 689 C If input structure hasn't been supplied from the PDB file read or generate
 690 C initial geometry.
 691         if (iranconf.eq.0 .and. .not. extconf) then
 692           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 693      &     write (iout,'(a)') 'Initial geometry will be read in.'
 694           if (read_cart) then
 695             read(inp,'(8f10.5)',end=36,err=36)
 696      &       ((c(l,k),l=1,3),k=1,nres),
 697      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 698             call int_from_cart1(.false.)
 699             do i=1,nres-1
 700               do j=1,3
 701                 dc(j,i)=c(j,i+1)-c(j,i)
 702                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 703               enddo
 704             enddo
 705             do i=nnt,nct
 706               if (itype(i).ne.10) then
 707                 do j=1,3
 708                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 709                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 710                 enddo
 711               endif
 712             enddo
 713             return
 714           else
 715             call read_angles(inp,*36)
 716           endif
 717           goto 37
 718    36     write (iout,'(a)') 'Error reading angle file.'
 719 #ifdef MPI
 720           call mpi_finalize( MPI_COMM_WORLD,IERR )
 721 #endif
 722           stop 'Error reading angle file.'
 723    37     continue
 724         else if (extconf) then
 725          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 726      &    write (iout,'(a)') 'Extended chain initial geometry.'
 727          do i=3,nres
 728           theta(i)=90d0*deg2rad
 729          enddo
 730          do i=4,nres
 731           phi(i)=180d0*deg2rad
 732          enddo
 733          do i=2,nres-1
 734           alph(i)=110d0*deg2rad
 735          enddo
 736          do i=2,nres-1
 737           omeg(i)=-120d0*deg2rad
 738          enddo
 739         else
 740           if(me.eq.king.or..not.out1file)
 741      &     write (iout,'(a)') 'Random-generated initial geometry.'
 742
 743
 744 #ifdef MPI
 745           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 746      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 747 #endif
 748             do itrial=1,100
 749               itmp=1
 750 c              call gen_rand_conf(itmp,*30)
 751               goto 40
 752    30         write (iout,*) 'Failed to generate random conformation',
 753      &          ', itrial=',itrial
 754               write (*,*) 'Processor:',me,
 755      &          ' Failed to generate random conformation',
 756      &          ' itrial=',itrial
 757               call intout
 758
 759 #ifdef AIX
 760               call flush_(iout)
 761 #else
 762               call flush(iout)
 763 #endif
 764             enddo
 765             write (iout,'(a,i3,a)') 'Processor:',me,
 766      &        ' error in generating random conformation.'
 767             write (*,'(a,i3,a)') 'Processor:',me,
 768      &        ' error in generating random conformation.'
 769             call flush(iout)
 770 #ifdef MPI
 771             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 772    40       continue
 773           endif
 774 #else
 775           do itrial=1,100
 776             itmp=1
 777 c            call gen_rand_conf(itmp,*31)
 778             goto 40
 779    31       write (iout,*) 'Failed to generate random conformation',
 780      &        ', itrial=',itrial
 781             write (*,*) 'Failed to generate random conformation',
 782      &        ', itrial=',itrial
 783           enddo
 784           write (iout,'(a,i3,a)') 'Processor:',me,
 785      &      ' error in generating random conformation.'
 786           write (*,'(a,i3,a)') 'Processor:',me,
 787      &      ' error in generating random conformation.'
 788           stop
 789    40     continue
 790 #endif
 791         endif
 792       elseif (modecalc.eq.4) then
 793         read (inp,'(a)') intinname
 794         open (intin,file=intinname,status='old',err=333)
 795         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 796      &  write (iout,'(a)') 'intinname',intinname
 797         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 798         goto 334
 799   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 800 #ifdef MPI
 801         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 802 #endif
 803         stop 'Error opening angle file.'
 804   334   continue
 805
 806       endif
 807 C Generate distance constraints, if the PDB structure is to be regularized.
 808       if (nthread.gt.0) then
 809         call read_threadbase
 810       endif
 811       call setup_var
 812       if (me.eq.king .or. .not. out1file)
 813      & call intout
 814       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 815         write (iout,'(/a,i3,a)')
 816      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 817         write (iout,'(20i4)') (iss(i),i=1,ns)
 818         write (iout,'(/a/)') 'Pre-formed links are:'
 819         do i=1,nss
 820           i1=ihpb(i)-nres
 821           i2=jhpb(i)-nres
 822           it1=itype(i1)
 823           it2=itype(i2)
 824           if (me.eq.king.or..not.out1file)
 825      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 826      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 827      &    ebr,forcon(i)
 828         enddo
 829         write (iout,'(a)')
 830       endif
 831 c       if (i2ndstr.gt.0) call secstrp2dihc
 832 c      call geom_to_var(nvar,x)
 833 c      call etotal(energia(0))
 834 c      call enerprint(energia(0))
 835 c      call briefout(0,etot)
 836 c      stop
 837 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 838 cd    write (iout,'(a)') 'Variable list:'
 839 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 840 #ifdef MPI
 841       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 842      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 843      &  'Processor',myrank,': end reading molecular data.'
 844 #endif
 845       return
 846       end
 847 c--------------------------------------------------------------------------
 848       logical function seq_comp(itypea,itypeb,length)
 849       implicit none
 850       integer length,itypea(length),itypeb(length)
 851       integer i
 852       do i=1,length
 853         if (itypea(i).ne.itypeb(i)) then
 854           seq_comp=.false.
 855           return
 856         endif
 857       enddo
 858       seq_comp=.true.
 859       return
 860       end
 861 c-----------------------------------------------------------------------------
 862       subroutine read_bridge
 863 C Read information about disulfide bridges.
 864       implicit real*8 (a-h,o-z)
 865       include 'DIMENSIONS'
 866 #ifdef MPI
 867       include 'mpif.h'
 868 #endif
 869       include 'COMMON.IOUNITS'
 870       include 'COMMON.GEO'
 871       include 'COMMON.VAR'
 872       include 'COMMON.INTERACT'
 873       include 'COMMON.LOCAL'
 874       include 'COMMON.NAMES'
 875       include 'COMMON.CHAIN'
 876       include 'COMMON.FFIELD'
 877       include 'COMMON.SBRIDGE'
 878       include 'COMMON.HEADER'
 879       include 'COMMON.CONTROL'
 880 c      include 'COMMON.DBASE'
 881 c      include 'COMMON.THREAD'
 882       include 'COMMON.TIME1'
 883       include 'COMMON.SETUP'
 884 C Read bridging residues.
 885       read (inp,*) ns,(iss(i),i=1,ns)
 886       print *,'ns=',ns
 887       if(me.eq.king.or..not.out1file)
 888      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 889 C Check whether the specified bridging residues are cystines.
 890       do i=1,ns
 891         if (itype(iss(i)).ne.1) then
 892           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 893      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 894      &   ' can form a disulfide bridge?!!!'
 895           write (*,'(2a,i3,a)')
 896      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 897      &   ' can form a disulfide bridge?!!!'
 898 #ifdef MPI
 899          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 900          stop
 901 #endif
 902         endif
 903       enddo
 904 C Read preformed bridges.
 905       if (ns.gt.0) then
 906       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 907       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 908       if (nss.gt.0) then
 909         nhpb=nss
 910 C Check if the residues involved in bridges are in the specified list of
 911 C bridging residues.
 912         do i=1,nss
 913           do j=1,i-1
 914             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 915      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 916               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 917      &      ' contains residues present in other pairs.'
 918               write (*,'(a,i3,a)') 'Disulfide pair',i,
 919      &      ' contains residues present in other pairs.'
 920 #ifdef MPI
 921               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 922               stop
 923 #endif
 924             endif
 925           enddo
 926           do j=1,ns
 927             if (ihpb(i).eq.iss(j)) goto 10
 928           enddo
 929           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 930    10     continue
 931           do j=1,ns
 932             if (jhpb(i).eq.iss(j)) goto 20
 933           enddo
 934           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 935    20     continue
 936           dhpb(i)=dbr
 937           forcon(i)=fbr
 938         enddo
 939         do i=1,nss
 940           ihpb(i)=ihpb(i)+nres
 941           jhpb(i)=jhpb(i)+nres
 942         enddo
 943       endif
 944       endif
 945       return
 946       end
 947 c----------------------------------------------------------------------------
 948       subroutine read_x(kanal,*)
 949       implicit real*8 (a-h,o-z)
 950       include 'DIMENSIONS'
 951       include 'COMMON.GEO'
 952       include 'COMMON.VAR'
 953       include 'COMMON.CHAIN'
 954       include 'COMMON.IOUNITS'
 955       include 'COMMON.CONTROL'
 956       include 'COMMON.LOCAL'
 957       include 'COMMON.INTERACT'
 958 c Read coordinates from input
 959 c
 960       read(kanal,'(8f10.5)',end=10,err=10)
 961      &  ((c(l,k),l=1,3),k=1,nres),
 962      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 963       do j=1,3
 964         c(j,nres+1)=c(j,1)
 965         c(j,2*nres)=c(j,nres)
 966       enddo
 967       call int_from_cart1(.false.)
 968       do i=1,nres-1
 969         do j=1,3
 970           dc(j,i)=c(j,i+1)-c(j,i)
 971           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 972         enddo
 973       enddo
 974       do i=nnt,nct
 975         if (itype(i).ne.10) then
 976           do j=1,3
 977             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 978             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 979           enddo
 980         endif
 981       enddo
 982
 983       return
 984    10 return1
 985       end
 986 c------------------------------------------------------------------------------
 987       subroutine setup_var
 988       implicit real*8 (a-h,o-z)
 989       include 'DIMENSIONS'
 990       include 'COMMON.IOUNITS'
 991       include 'COMMON.GEO'
 992       include 'COMMON.VAR'
 993       include 'COMMON.INTERACT'
 994       include 'COMMON.LOCAL'
 995       include 'COMMON.NAMES'
 996       include 'COMMON.CHAIN'
 997       include 'COMMON.FFIELD'
 998       include 'COMMON.SBRIDGE'
 999       include 'COMMON.HEADER'
1000       include 'COMMON.CONTROL'
1001 c      include 'COMMON.DBASE'
1002 c      include 'COMMON.THREAD'
1003       include 'COMMON.TIME1'
1004 C Set up variable list.
1005       ntheta=nres-2
1006       nphi=nres-3
1007       nvar=ntheta+nphi
1008       nside=0
1009       do i=2,nres-1
1010         if (itype(i).ne.10) then
1011           nside=nside+1
1012           ialph(i,1)=nvar+nside
1013           ialph(nside,2)=i
1014         endif
1015       enddo
1016       if (indphi.gt.0) then
1017         nvar=nphi
1018       else if (indback.gt.0) then
1019         nvar=nphi+ntheta
1020       else
1021         nvar=nvar+2*nside
1022       endif
1023 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1024       return
1025       end
1026 c----------------------------------------------------------------------------
1027       subroutine gen_dist_constr
1028 C Generate CA distance constraints.
1029       implicit real*8 (a-h,o-z)
1030       include 'DIMENSIONS'
1031       include 'COMMON.IOUNITS'
1032       include 'COMMON.GEO'
1033       include 'COMMON.VAR'
1034       include 'COMMON.INTERACT'
1035       include 'COMMON.LOCAL'
1036       include 'COMMON.NAMES'
1037       include 'COMMON.CHAIN'
1038       include 'COMMON.FFIELD'
1039       include 'COMMON.SBRIDGE'
1040       include 'COMMON.HEADER'
1041       include 'COMMON.CONTROL'
1042 c      include 'COMMON.DBASE'
1043 c      include 'COMMON.THREAD'
1044       include 'COMMON.TIME1'
1045       dimension itype_pdb(maxres)
1046       common /pizda/ itype_pdb
1047       character*2 iden
1048 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1049 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1050 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1051 cd     & ' nsup',nsup
1052       do i=nstart_sup,nstart_sup+nsup-1
1053 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1054 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1055         do j=i+2,nstart_sup+nsup-1
1056           nhpb=nhpb+1
1057           ihpb(nhpb)=i+nstart_seq-nstart_sup
1058           jhpb(nhpb)=j+nstart_seq-nstart_sup
1059           forcon(nhpb)=weidis
1060           dhpb(nhpb)=dist(i,j)
1061         enddo
1062       enddo
1063 cd      write (iout,'(a)') 'Distance constraints:'
1064 cd      do i=nss+1,nhpb
1065 cd        ii=ihpb(i)
1066 cd        jj=jhpb(i)
1067 cd        iden='CA'
1068 cd        if (ii.gt.nres) then
1069 cd          iden='SC'
1070 cd          ii=ii-nres
1071 cd          jj=jj-nres
1072 cd        endif
1073 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1074 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1075 cd     &  dhpb(i),forcon(i)
1076 cd      enddo
1077       return
1078       end
1079 c----------------------------------------------------------------------------
1080       subroutine csaread
1081       implicit real*8 (a-h,o-z)
1082       include 'DIMENSIONS'
1083       include 'COMMON.IOUNITS'
1084       include 'COMMON.GEO'
1085       include 'COMMON.CSA'
1086       include 'COMMON.BANK'
1087       include 'COMMON.CONTROL'
1088       character*80 ucase
1089       character*620 mcmcard
1090       call card_concat(mcmcard)
1091
1092       call readi(mcmcard,'NCONF',nconf,50)
1093       call readi(mcmcard,'NADD',nadd,0)
1094       call readi(mcmcard,'JSTART',jstart,1)
1095       call readi(mcmcard,'JEND',jend,1)
1096       call readi(mcmcard,'NSTMAX',nstmax,500000)
1097       call readi(mcmcard,'N0',n0,1)
1098       call readi(mcmcard,'N1',n1,6)
1099       call readi(mcmcard,'N2',n2,4)
1100       call readi(mcmcard,'N3',n3,0)
1101       call readi(mcmcard,'N4',n4,0)
1102       call readi(mcmcard,'N5',n5,0)
1103       call readi(mcmcard,'N6',n6,10)
1104       call readi(mcmcard,'N7',n7,0)
1105       call readi(mcmcard,'N8',n8,0)
1106       call readi(mcmcard,'N9',n9,0)
1107       call readi(mcmcard,'N14',n14,0)
1108       call readi(mcmcard,'N15',n15,0)
1109       call readi(mcmcard,'N16',n16,0)
1110       call readi(mcmcard,'N17',n17,0)
1111       call readi(mcmcard,'N18',n18,0)
1112
1113       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1114
1115       call readi(mcmcard,'NDIFF',ndiff,2)
1116       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1117       call readi(mcmcard,'IS1',is1,1)
1118       call readi(mcmcard,'IS2',is2,8)
1119       call readi(mcmcard,'NRAN0',nran0,4)
1120       call readi(mcmcard,'NRAN1',nran1,2)
1121       call readi(mcmcard,'IRR',irr,1)
1122       call readi(mcmcard,'NSEED',nseed,20)
1123       call readi(mcmcard,'NTOTAL',ntotal,10000)
1124       call reada(mcmcard,'CUT1',cut1,2.0d0)
1125       call reada(mcmcard,'CUT2',cut2,5.0d0)
1126       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1127       call readi(mcmcard,'ICMAX',icmax,1)
1128       call readi(mcmcard,'NBANKM',nbankm,400)
1129       call readi(mcmcard,'IUCUT',iucut,2)
1130       call readi(mcmcard,'IRESTART',irestart,0)
1131 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1132       ntbankm=0
1133 c!bankt
1134       call reada(mcmcard,'DELE',dele,20.0d0)
1135       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1136       call readi(mcmcard,'IREF',iref,0)
1137       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1138       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1139       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1140       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1141       write (iout,*) "NCONF_IN",nconf_in
1142       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1143       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1144      &      " for torsional angles"
1145       return
1146       end
1147
1148       subroutine read_minim
1149       implicit real*8 (a-h,o-z)
1150       include 'DIMENSIONS'
1151       include 'COMMON.MINIM'
1152       include 'COMMON.IOUNITS'
1153       character*80 ucase
1154       character*320 minimcard
1155       call card_concat(minimcard)
1156       call readi(minimcard,'MAXMIN',maxmin,2000)
1157       call readi(minimcard,'MAXFUN',maxfun,5000)
1158       call readi(minimcard,'MINMIN',minmin,maxmin)
1159       call readi(minimcard,'MINFUN',minfun,maxmin)
1160       call reada(minimcard,'TOLF',tolf,1.0D-2)
1161       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1162       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1163       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1164       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1165       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1166      &         'Options in energy minimization:'
1167       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1168      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1169      & 'MinMin:',MinMin,' MinFun:',MinFun,
1170      & ' TolF:',TolF,' RTolF:',RTolF
1171       return
1172       end
1173 c----------------------------------------------------------------------------
1174       subroutine read_angles(kanal,*)
1175       implicit real*8 (a-h,o-z)
1176       include 'DIMENSIONS'
1177       include 'COMMON.GEO'
1178       include 'COMMON.VAR'
1179       include 'COMMON.CHAIN'
1180       include 'COMMON.IOUNITS'
1181       include 'COMMON.CONTROL'
1182 c Read angles from input
1183 c
1184        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1185        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1186        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1187        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1188
1189        do i=1,nres
1190 c 9/7/01 avoid 180 deg valence angle
1191         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1192 c
1193         theta(i)=deg2rad*theta(i)
1194         phi(i)=deg2rad*phi(i)
1195         alph(i)=deg2rad*alph(i)
1196         omeg(i)=deg2rad*omeg(i)
1197        enddo
1198       return
1199    10 return1
1200       end
1201 c----------------------------------------------------------------------------
1202       subroutine reada(rekord,lancuch,wartosc,default)
1203       implicit none
1204       character*(*) rekord,lancuch
1205       double precision wartosc,default
1206       integer ilen,iread
1207       external ilen
1208       iread=index(rekord,lancuch)
1209       if (iread.eq.0) then
1210         wartosc=default
1211         return
1212       endif
1213       iread=iread+ilen(lancuch)+1
1214       read (rekord(iread:),*,err=10,end=10) wartosc
1215       return
1216   10  wartosc=default
1217       return
1218       end
1219 c----------------------------------------------------------------------------
1220       subroutine readi(rekord,lancuch,wartosc,default)
1221       implicit none
1222       character*(*) rekord,lancuch
1223       integer wartosc,default
1224       integer ilen,iread
1225       external ilen
1226       iread=index(rekord,lancuch)
1227       if (iread.eq.0) then
1228         wartosc=default
1229         return
1230       endif
1231       iread=iread+ilen(lancuch)+1
1232       read (rekord(iread:),*,err=10,end=10) wartosc
1233       return
1234   10  wartosc=default
1235       return
1236       end
1237 c----------------------------------------------------------------------------
1238       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1239       implicit none
1240       integer dim,i
1241       integer tablica(dim),default
1242       character*(*) rekord,lancuch
1243       character*80 aux
1244       integer ilen,iread
1245       external ilen
1246       do i=1,dim
1247         tablica(i)=default
1248       enddo
1249       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1250       if (iread.eq.0) return
1251       iread=iread+ilen(lancuch)+1
1252       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1253    10 return
1254       end
1255 c----------------------------------------------------------------------------
1256       subroutine multreada(rekord,lancuch,tablica,dim,default)
1257       implicit none
1258       integer dim,i
1259       double precision tablica(dim),default
1260       character*(*) rekord,lancuch
1261       character*80 aux
1262       integer ilen,iread
1263       external ilen
1264       do i=1,dim
1265         tablica(i)=default
1266       enddo
1267       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1268       if (iread.eq.0) return
1269       iread=iread+ilen(lancuch)+1
1270       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1271    10 return
1272       end
1273 c----------------------------------------------------------------------------
1274       subroutine openunits
1275       implicit real*8 (a-h,o-z)
1276       include 'DIMENSIONS'
1277 #ifdef MPI
1278       include 'mpif.h'
1279       character*16 form,nodename
1280       integer nodelen
1281 #endif
1282       include 'COMMON.SETUP'
1283       include 'COMMON.IOUNITS'
1284 c      include 'COMMON.MD'
1285       include 'COMMON.CONTROL'
1286       integer lenpre,lenpot,ilen,lentmp
1287       external ilen
1288       character*3 out1file_text,ucase
1289       character*3 ll
1290       external ucase
1291 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1292       call getenv_loc("PREFIX",prefix)
1293       pref_orig = prefix
1294       call getenv_loc("POT",pot)
1295       call getenv_loc("DIRTMP",tmpdir)
1296       call getenv_loc("CURDIR",curdir)
1297       call getenv_loc("OUT1FILE",out1file_text)
1298 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1299       out1file_text=ucase(out1file_text)
1300       if (out1file_text(1:1).eq."Y") then
1301         out1file=.true.
1302       else
1303         out1file=fg_rank.gt.0
1304       endif
1305       lenpre=ilen(prefix)
1306       lenpot=ilen(pot)
1307       lentmp=ilen(tmpdir)
1308       if (lentmp.gt.0) then
1309           write (*,'(80(1h!))')
1310           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1311           write (*,'(80(1h!))')
1312           write (*,*)"All output files will be on node /tmp directory."
1313 #ifdef MPI
1314         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1315         if (me.eq.king) then
1316           write (*,*) "The master node is ",nodename
1317         else if (fg_rank.eq.0) then
1318           write (*,*) "I am the CG slave node ",nodename
1319         else
1320           write (*,*) "I am the FG slave node ",nodename
1321         endif
1322 #endif
1323         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1324         lenpre = lentmp+lenpre+1
1325       endif
1326       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1327 C Get the names and open the input files
1328 #if defined(WINIFL) || defined(WINPGI)
1329       open(1,file=pref_orig(:ilen(pref_orig))//
1330      &  '.inp',status='old',readonly,shared)
1331        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1332 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1333 C Get parameter filenames and open the parameter files.
1334       call getenv_loc('BONDPAR',bondname)
1335       open (ibond,file=bondname,status='old',readonly,shared)
1336       call getenv_loc('THETPAR',thetname)
1337       open (ithep,file=thetname,status='old',readonly,shared)
1338 #ifndef CRYST_THETA
1339       call getenv_loc('THETPARPDB',thetname_pdb)
1340       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1341 #endif
1342       call getenv_loc('ROTPAR',rotname)
1343       open (irotam,file=rotname,status='old',readonly,shared)
1344 #ifndef CRYST_SC
1345       call getenv_loc('ROTPARPDB',rotname_pdb)
1346       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1347 #endif
1348       call getenv_loc('TORPAR',torname)
1349       open (itorp,file=torname,status='old',readonly,shared)
1350       call getenv_loc('TORDPAR',tordname)
1351       open (itordp,file=tordname,status='old',readonly,shared)
1352       call getenv_loc('FOURIER',fouriername)
1353       open (ifourier,file=fouriername,status='old',readonly,shared)
1354       call getenv_loc('ELEPAR',elename)
1355       open (ielep,file=elename,status='old',readonly,shared)
1356       call getenv_loc('SIDEPAR',sidename)
1357       open (isidep,file=sidename,status='old',readonly,shared)
1358 #elif (defined CRAY) || (defined AIX)
1359       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1360      &  action='read')
1361 c      print *,"Processor",myrank," opened file 1"
1362       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1363 c      print *,"Processor",myrank," opened file 9"
1364 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1365 C Get parameter filenames and open the parameter files.
1366       call getenv_loc('BONDPAR',bondname)
1367       open (ibond,file=bondname,status='old',action='read')
1368 c      print *,"Processor",myrank," opened file IBOND"
1369       call getenv_loc('THETPAR',thetname)
1370       open (ithep,file=thetname,status='old',action='read')
1371 c      print *,"Processor",myrank," opened file ITHEP"
1372 #ifndef CRYST_THETA
1373       call getenv_loc('THETPARPDB',thetname_pdb)
1374       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1375 #endif
1376       call getenv_loc('ROTPAR',rotname)
1377       open (irotam,file=rotname,status='old',action='read')
1378 c      print *,"Processor",myrank," opened file IROTAM"
1379 #ifndef CRYST_SC
1380       call getenv_loc('ROTPARPDB',rotname_pdb)
1381       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1382 #endif
1383       call getenv_loc('TORPAR',torname)
1384       open (itorp,file=torname,status='old',action='read')
1385 c      print *,"Processor",myrank," opened file ITORP"
1386       call getenv_loc('TORDPAR',tordname)
1387       open (itordp,file=tordname,status='old',action='read')
1388 c      print *,"Processor",myrank," opened file ITORDP"
1389       call getenv_loc('SCCORPAR',sccorname)
1390       open (isccor,file=sccorname,status='old',action='read')
1391 c      print *,"Processor",myrank," opened file ISCCOR"
1392       call getenv_loc('FOURIER',fouriername)
1393       open (ifourier,file=fouriername,status='old',action='read')
1394 c      print *,"Processor",myrank," opened file IFOURIER"
1395       call getenv_loc('ELEPAR',elename)
1396       open (ielep,file=elename,status='old',action='read')
1397 c      print *,"Processor",myrank," opened file IELEP"
1398       call getenv_loc('SIDEPAR',sidename)
1399       open (isidep,file=sidename,status='old',action='read')
1400 c      print *,"Processor",myrank," opened file ISIDEP"
1401 c      print *,"Processor",myrank," opened parameter files"
1402 #elif (defined G77)
1403       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1404       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1405 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1406 C Get parameter filenames and open the parameter files.
1407       call getenv_loc('BONDPAR',bondname)
1408       open (ibond,file=bondname,status='old')
1409       call getenv_loc('THETPAR',thetname)
1410       open (ithep,file=thetname,status='old')
1411 #ifndef CRYST_THETA
1412       call getenv_loc('THETPARPDB',thetname_pdb)
1413       open (ithep_pdb,file=thetname_pdb,status='old')
1414 #endif
1415       call getenv_loc('ROTPAR',rotname)
1416       open (irotam,file=rotname,status='old')
1417 #ifndef CRYST_SC
1418       call getenv_loc('ROTPARPDB',rotname_pdb)
1419       open (irotam_pdb,file=rotname_pdb,status='old')
1420 #endif
1421       call getenv_loc('TORPAR',torname)
1422       open (itorp,file=torname,status='old')
1423       call getenv_loc('TORDPAR',tordname)
1424       open (itordp,file=tordname,status='old')
1425       call getenv_loc('SCCORPAR',sccorname)
1426       open (isccor,file=sccorname,status='old')
1427       call getenv_loc('FOURIER',fouriername)
1428       open (ifourier,file=fouriername,status='old')
1429       call getenv_loc('ELEPAR',elename)
1430       open (ielep,file=elename,status='old')
1431       call getenv_loc('SIDEPAR',sidename)
1432       open (isidep,file=sidename,status='old')
1433 #else
1434       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1435      &  readonly)
1436        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1437 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1438 C Get parameter filenames and open the parameter files.
1439       call getenv_loc('BONDPAR',bondname)
1440       open (ibond,file=bondname,status='old',readonly)
1441       call getenv_loc('THETPAR',thetname)
1442       open (ithep,file=thetname,status='old',readonly)
1443 #ifndef CRYST_THETA
1444       call getenv_loc('THETPARPDB',thetname_pdb)
1445       print *,"thetname_pdb ",thetname_pdb
1446       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1447       print *,ithep_pdb," opened"
1448 #endif
1449       call getenv_loc('ROTPAR',rotname)
1450       open (irotam,file=rotname,status='old',readonly)
1451 #ifndef CRYST_SC
1452       call getenv_loc('ROTPARPDB',rotname_pdb)
1453       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1454 #endif
1455       call getenv_loc('TORPAR',torname)
1456       open (itorp,file=torname,status='old',readonly)
1457       call getenv_loc('TORDPAR',tordname)
1458       open (itordp,file=tordname,status='old',readonly)
1459       call getenv_loc('SCCORPAR',sccorname)
1460       open (isccor,file=sccorname,status='old',readonly)
1461       call getenv_loc('FOURIER',fouriername)
1462       open (ifourier,file=fouriername,status='old',readonly)
1463       call getenv_loc('ELEPAR',elename)
1464       open (ielep,file=elename,status='old',readonly)
1465       call getenv_loc('SIDEPAR',sidename)
1466       open (isidep,file=sidename,status='old',readonly)
1467 #endif
1468 #ifndef OLDSCP
1469 C
1470 C 8/9/01 In the newest version SCp interaction constants are read from a file
1471 C Use -DOLDSCP to use hard-coded constants instead.
1472 C
1473       call getenv_loc('SCPPAR',scpname)
1474 #if defined(WINIFL) || defined(WINPGI)
1475       open (iscpp,file=scpname,status='old',readonly,shared)
1476 #elif (defined CRAY)  || (defined AIX)
1477       open (iscpp,file=scpname,status='old',action='read')
1478 #elif (defined G77)
1479       open (iscpp,file=scpname,status='old')
1480 #else
1481       open (iscpp,file=scpname,status='old',readonly)
1482 #endif
1483 #endif
1484       call getenv_loc('PATTERN',patname)
1485 #if defined(WINIFL) || defined(WINPGI)
1486       open (icbase,file=patname,status='old',readonly,shared)
1487 #elif (defined CRAY)  || (defined AIX)
1488       open (icbase,file=patname,status='old',action='read')
1489 #elif (defined G77)
1490       open (icbase,file=patname,status='old')
1491 #else
1492       open (icbase,file=patname,status='old',readonly)
1493 #endif
1494 #ifdef MPI
1495 C Open output file only for CG processes
1496 c      print *,"Processor",myrank," fg_rank",fg_rank
1497       if (fg_rank.eq.0) then
1498
1499       if (nodes.eq.1) then
1500         npos=3
1501       else
1502         npos = dlog10(dfloat(nodes-1))+1
1503       endif
1504       if (npos.lt.3) npos=3
1505       write (liczba,'(i1)') npos
1506       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1507      &  //')'
1508       write (liczba,form) me
1509       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1510      &  liczba(:ilen(liczba))
1511       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1512      &  //'.int'
1513       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1514      &  //'.pdb'
1515       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1516      &  liczba(:ilen(liczba))//'.mol2'
1517       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1518      &  liczba(:ilen(liczba))//'.stat'
1519       if (lentmp.gt.0)
1520      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1521      &      //liczba(:ilen(liczba))//'.stat')
1522       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1523      &  //'.rst'
1524 c      if(usampl) then
1525 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1526 c     & liczba(:ilen(liczba))//'.const'
1527 c      endif
1528
1529       endif
1530 #else
1531       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1532       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1533       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1534       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1535       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1536       if (lentmp.gt.0)
1537      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1538      &    '.stat')
1539       rest2name=prefix(:ilen(prefix))//'.rst'
1540 c      if(usampl) then
1541 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1542 c      endif
1543 #endif
1544 #if defined(AIX) || defined(PGI)
1545       if (me.eq.king .or. .not. out1file)
1546      &   open(iout,file=outname,status='unknown')
1547 #ifdef DEBUG
1548       if (fg_rank.gt.0) then
1549         write (liczba,'(i3.3)') myrank/nfgtasks
1550         write (ll,'(bz,i3.3)') fg_rank
1551         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1552      &   status='unknown')
1553       endif
1554 #endif
1555       if(me.eq.king) then
1556        open(igeom,file=intname,status='unknown',position='append')
1557        open(ipdb,file=pdbname,status='unknown')
1558        open(imol2,file=mol2name,status='unknown')
1559        open(istat,file=statname,status='unknown',position='append')
1560       else
1561 c1out       open(iout,file=outname,status='unknown')
1562       endif
1563 #else
1564       if (me.eq.king .or. .not.out1file)
1565      &    open(iout,file=outname,status='unknown')
1566 #ifdef DEBUG
1567       if (fg_rank.gt.0) then
1568         write (liczba,'(i3.3)') myrank/nfgtasks
1569         write (ll,'(bz,i3.3)') fg_rank
1570         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1571      &   status='unknown')
1572       endif
1573 #endif
1574       if(me.eq.king) then
1575        open(igeom,file=intname,status='unknown',access='append')
1576        open(ipdb,file=pdbname,status='unknown')
1577        open(imol2,file=mol2name,status='unknown')
1578        open(istat,file=statname,status='unknown',access='append')
1579       else
1580 c1out       open(iout,file=outname,status='unknown')
1581       endif
1582 #endif
1583       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1584       csa_seed=prefix(:lenpre)//'.CSA.seed'
1585       csa_history=prefix(:lenpre)//'.CSA.history'
1586       csa_bank=prefix(:lenpre)//'.CSA.bank'
1587       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1588       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1589       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1590 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1591       csa_int=prefix(:lenpre)//'.int'
1592       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1593       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1594       csa_in=prefix(:lenpre)//'.CSA.in'
1595 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1596 C Write file names
1597       if (me.eq.king)then
1598       write (iout,'(80(1h-))')
1599       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1600       write (iout,'(80(1h-))')
1601       write (iout,*) "Input file                      : ",
1602      &  pref_orig(:ilen(pref_orig))//'.inp'
1603       write (iout,*) "Output file                     : ",
1604      &  outname(:ilen(outname))
1605       write (iout,*)
1606       write (iout,*) "Sidechain potential file        : ",
1607      &  sidename(:ilen(sidename))
1608 #ifndef OLDSCP
1609       write (iout,*) "SCp potential file              : ",
1610      &  scpname(:ilen(scpname))
1611 #endif
1612       write (iout,*) "Electrostatic potential file    : ",
1613      &  elename(:ilen(elename))
1614       write (iout,*) "Cumulant coefficient file       : ",
1615      &  fouriername(:ilen(fouriername))
1616       write (iout,*) "Torsional parameter file        : ",
1617      &  torname(:ilen(torname))
1618       write (iout,*) "Double torsional parameter file : ",
1619      &  tordname(:ilen(tordname))
1620       write (iout,*) "SCCOR parameter file : ",
1621      &  sccorname(:ilen(sccorname))
1622       write (iout,*) "Bond & inertia constant file    : ",
1623      &  bondname(:ilen(bondname))
1624       write (iout,*) "Bending parameter file          : ",
1625      &  thetname(:ilen(thetname))
1626       write (iout,*) "Rotamer parameter file          : ",
1627      &  rotname(:ilen(rotname))
1628       write (iout,*) "Threading database              : ",
1629      &  patname(:ilen(patname))
1630       if (lentmp.ne.0)
1631      &write (iout,*)" DIRTMP                          : ",
1632      &  tmpdir(:lentmp)
1633       write (iout,'(80(1h-))')
1634       endif
1635       return
1636       end
1637 c----------------------------------------------------------------------------
1638       subroutine card_concat(card)
1639       implicit real*8 (a-h,o-z)
1640       include 'DIMENSIONS'
1641       include 'COMMON.IOUNITS'
1642       character*(*) card
1643       character*80 karta,ucase
1644       external ilen
1645       read (inp,'(a)') karta
1646       karta=ucase(karta)
1647       card=' '
1648       do while (karta(80:80).eq.'&')
1649         card=card(:ilen(card)+1)//karta(:79)
1650         read (inp,'(a)') karta
1651         karta=ucase(karta)
1652       enddo
1653       card=card(:ilen(card)+1)//karta
1654       return
1655       end
1656
1657       subroutine read_dist_constr
1658       implicit real*8 (a-h,o-z)
1659       include 'DIMENSIONS'
1660 #ifdef MPI
1661       include 'mpif.h'
1662 #endif
1663       include 'COMMON.SETUP'
1664       include 'COMMON.CONTROL'
1665       include 'COMMON.CHAIN'
1666       include 'COMMON.IOUNITS'
1667       include 'COMMON.SBRIDGE'
1668       integer ifrag_(2,100),ipair_(2,100)
1669       double precision wfrag_(100),wpair_(100)
1670       character*500 controlcard
1671       write (iout,*) "Calling read_dist_constr"
1672       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1673       call flush(iout)
1674       call card_concat(controlcard)
1675       call readi(controlcard,"NFRAG",nfrag_,0)
1676       call readi(controlcard,"NPAIR",npair_,0)
1677       call readi(controlcard,"NDIST",ndist_,0)
1678       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1679       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1680       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1681       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1682       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1683 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1684 c      write (iout,*) "IFRAG"
1685 c      do i=1,nfrag_
1686 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1687 c      enddo
1688 c      write (iout,*) "IPAIR"
1689 c      do i=1,npair_
1690 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1691 c      enddo
1692       call flush(iout)
1693       do i=1,nfrag_
1694         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1695         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1696      &    ifrag_(2,i)=nstart_sup+nsup-1
1697 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1698         call flush(iout)
1699         if (wfrag_(i).gt.0.0d0) then
1700         do j=ifrag_(1,i),ifrag_(2,i)-1
1701           do k=j+1,ifrag_(2,i)
1702             write (iout,*) "j",j," k",k
1703             ddjk=dist(j,k)
1704             if (constr_dist.eq.1) then
1705             nhpb=nhpb+1
1706             ihpb(nhpb)=j
1707             jhpb(nhpb)=k
1708               dhpb(nhpb)=ddjk
1709             forcon(nhpb)=wfrag_(i)
1710             else if (constr_dist.eq.2) then
1711               if (ddjk.le.dist_cut) then
1712                 nhpb=nhpb+1
1713                 ihpb(nhpb)=j
1714                 jhpb(nhpb)=k
1715                 dhpb(nhpb)=ddjk
1716                 forcon(nhpb)=wfrag_(i)
1717               endif
1718             else
1719               nhpb=nhpb+1
1720               ihpb(nhpb)=j
1721               jhpb(nhpb)=k
1722               dhpb(nhpb)=ddjk
1723               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1724             endif
1725 #ifdef MPI
1726             if (.not.out1file .or. me.eq.king)
1727      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1728      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1729 #else
1730             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1731      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1732 #endif
1733           enddo
1734         enddo
1735         endif
1736       enddo
1737       do i=1,npair_
1738         if (wpair_(i).gt.0.0d0) then
1739         ii = ipair_(1,i)
1740         jj = ipair_(2,i)
1741         if (ii.gt.jj) then
1742           itemp=ii
1743           ii=jj
1744           jj=itemp
1745         endif
1746         do j=ifrag_(1,ii),ifrag_(2,ii)
1747           do k=ifrag_(1,jj),ifrag_(2,jj)
1748             nhpb=nhpb+1
1749             ihpb(nhpb)=j
1750             jhpb(nhpb)=k
1751             forcon(nhpb)=wpair_(i)
1752             dhpb(nhpb)=dist(j,k)
1753 #ifdef MPI
1754             if (.not.out1file .or. me.eq.king)
1755      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1756      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1757 #else
1758             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1759      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1760 #endif
1761           enddo
1762         enddo
1763         endif
1764       enddo
1765       do i=1,ndist_
1766         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1767         if (forcon(nhpb+1).gt.0.0d0) then
1768           nhpb=nhpb+1
1769           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1770 #ifdef MPI
1771           if (.not.out1file .or. me.eq.king)
1772      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1773      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1774 #else
1775           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1776      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1777 #endif
1778         endif
1779       enddo
1780       call flush(iout)
1781       return
1782       end
1783 c-------------------------------------------------------------------------------
1784 #ifdef WINIFL
1785       subroutine flush(iu)
1786       return
1787       end
1788 #endif
1789 #ifdef AIX
1790       subroutine flush(iu)
1791       call flush_(iu)
1792       return
1793       end
1794 #endif
1795 c------------------------------------------------------------------------------
1796       subroutine copy_to_tmp(source)
1797       include "DIMENSIONS"
1798       include "COMMON.IOUNITS"
1799       character*(*) source
1800       character* 256 tmpfile
1801       integer ilen
1802       external ilen
1803       logical ex
1804       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1805       inquire(file=tmpfile,exist=ex)
1806       if (ex) then
1807         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1808      &   " to temporary directory..."
1809         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1810         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1811       endif
1812       return
1813       end
1814 c------------------------------------------------------------------------------
1815       subroutine move_from_tmp(source)
1816       include "DIMENSIONS"
1817       include "COMMON.IOUNITS"
1818       character*(*) source
1819       integer ilen
1820       external ilen
1821       write (*,*) "Moving ",source(:ilen(source)),
1822      & " from temporary directory to working directory"
1823       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1824       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1825       return
1826       end
1827 c------------------------------------------------------------------------------
1828       subroutine random_init(seed)
1829 C
1830 C Initialize random number generator
1831 C
1832       implicit real*8 (a-h,o-z)
1833       include 'DIMENSIONS'
1834 #ifdef AMD64
1835       integer*8 iseedi8
1836 #endif
1837 #ifdef MPI
1838       include 'mpif.h'
1839       logical OKRandom, prng_restart
1840       real*8  r1
1841       integer iseed_array(4)
1842 #endif
1843       include 'COMMON.IOUNITS'
1844       include 'COMMON.TIME1'
1845 c      include 'COMMON.THREAD'
1846       include 'COMMON.SBRIDGE'
1847       include 'COMMON.CONTROL'
1848       include 'COMMON.MCM'
1849 c      include 'COMMON.MAP'
1850       include 'COMMON.HEADER'
1851 c      include 'COMMON.CSA'
1852       include 'COMMON.CHAIN'
1853 c      include 'COMMON.MUCA'
1854 c      include 'COMMON.MD'
1855       include 'COMMON.FFIELD'
1856       include 'COMMON.SETUP'
1857       iseed=-dint(dabs(seed))
1858       if (iseed.eq.0) then
1859         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1860      &    'Random seed undefined. The program will stop.'
1861         write (*,'(/80(1h*)/20x,a/80(1h*))')
1862      &    'Random seed undefined. The program will stop.'
1863 #ifdef MPI
1864         call mpi_finalize(mpi_comm_world,ierr)
1865 #endif
1866         stop 'Bad random seed.'
1867       endif
1868 #ifdef MPI
1869       if (fg_rank.eq.0) then
1870       seed=seed*(me+1)+1
1871 #ifdef AMD64
1872       iseedi8=dint(seed)
1873       if(me.eq.king .or. .not. out1file)
1874      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1875       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1876       OKRandom = prng_restart(me,iseedi8)
1877 #else
1878       do i=1,4
1879        tmp=65536.0d0**(4-i)
1880        iseed_array(i) = dint(seed/tmp)
1881        seed=seed-iseed_array(i)*tmp
1882       enddo
1883       if(me.eq.king .or. .not. out1file)
1884      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1885      &                 (iseed_array(i),i=1,4)
1886       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1887      &                 (iseed_array(i),i=1,4)
1888       OKRandom = prng_restart(me,iseed_array)
1889 #endif
1890       if (OKRandom) then
1891         r1=ran_number(0.0D0,1.0D0)
1892         if(me.eq.king .or. .not. out1file)
1893      &   write (iout,*) 'ran_num',r1
1894         if (r1.lt.0.0d0) OKRandom=.false.
1895       endif
1896       if (.not.OKRandom) then
1897         write (iout,*) 'PRNG IS NOT WORKING!!!'
1898         print *,'PRNG IS NOT WORKING!!!'
1899         if (me.eq.0) then
1900          call flush(iout)
1901          call mpi_abort(mpi_comm_world,error_msg,ierr)
1902          stop
1903         else
1904          write (iout,*) 'too many processors for parallel prng'
1905          write (*,*) 'too many processors for parallel prng'
1906          call flush(iout)
1907          stop
1908         endif
1909       endif
1910       endif
1911 #else
1912       call vrndst(iseed)
1913       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1914 #endif
1915       return
1916       end