source/unres/src_CSA_DiL/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 C Following 2 lines for diagnostics; comment out if not needed
 459          write (iout,*) "Before sideadd"
 460          call intout
 461          if(me.eq.king.or..not.out1file)
 462      &    write(iout,*)'Adding sidechains'
 463          maxsi=1000
 464          do i=2,nres-1
 465           iti=itype(i)
 466           if (iti.ne.10) then
 467             nsi=0
 468             fail=.true.
 469             do while (fail.and.nsi.le.maxsi)
 470 c              call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 471               nsi=nsi+1
 472             enddo
 473             if(fail) write(iout,*)'Adding sidechain failed for res ',
 474      &              i,' after ',nsi,' trials'
 475           endif
 476          enddo
 477         endif
 478 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
 479         call chainbuild
 480 C Following 2 lines for diagnostics; comment out if not needed
 481         write (iout,*) "After sideadd"
 482         call intout
 483       endif
 484
 485       if (indpdb.eq.0) then
 486 C Read sequence if not taken from the pdb file.
 487         read (inp,*) nres
 488 c        print *,'nres=',nres
 489         if (iscode.gt.0) then
 490           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 491         else
 492           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 493         endif
 494 C Convert sequence to numeric code
 495         do i=1,nres
 496           itype(i)=rescode(i,sequence(i),iscode)
 497         enddo
 498 C Assign initial virtual bond lengths
 499         do i=2,nres
 500           vbld(i)=vbl
 501           vbld_inv(i)=vblinv
 502         enddo
 503         do i=2,nres-1
 504           vbld(i+nres)=dsc(iabs(itype(i)))
 505           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 506 c          write (iout,*) "i",i," itype",itype(i),
 507 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 508         enddo
 509       endif
 510 c      print *,nres
 511 c      print '(20i4)',(itype(i),i=1,nres)
 512       do i=1,nres
 513 #ifdef PROCOR
 514         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 515 #else
 516         if (itype(i).eq.21) then
 517 #endif
 518           itel(i)=0
 519 #ifdef PROCOR
 520         else if (iabs(itype(i+1)).ne.20) then
 521 #else
 522         else if (iabs(itype(i)).ne.20) then
 523 #endif
 524           itel(i)=1
 525         else
 526           itel(i)=2
 527         endif
 528       enddo
 529       if(me.eq.king.or..not.out1file)then
 530        write (iout,*) "ITEL"
 531        do i=1,nres-1
 532          write (iout,*) i,itype(i),itel(i)
 533        enddo
 534        print *,'Call Read_Bridge.'
 535       endif
 536       call read_bridge
 537 C 8/13/98 Set limits to generating the dihedral angles
 538       do i=1,nres
 539         phibound(1,i)=-pi
 540         phibound(2,i)=pi
 541       enddo
 542       read (inp,*) ndih_constr
 543       if (ndih_constr.gt.0) then
 544         read (inp,*) ftors
 545         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 546         if(me.eq.king.or..not.out1file)then
 547          write (iout,*)
 548      &   'There are',ndih_constr,' constraints on phi angles.'
 549          do i=1,ndih_constr
 550           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 551          enddo
 552         endif
 553         do i=1,ndih_constr
 554           phi0(i)=deg2rad*phi0(i)
 555           drange(i)=deg2rad*drange(i)
 556         enddo
 557         if(me.eq.king.or..not.out1file)
 558      &   write (iout,*) 'FTORS',ftors
 559         do i=1,ndih_constr
 560           ii = idih_constr(i)
 561           phibound(1,ii) = phi0(i)-drange(i)
 562           phibound(2,ii) = phi0(i)+drange(i)
 563         enddo
 564       endif
 565       nnt=1
 566 #ifdef MPI
 567       if (me.eq.king) then
 568 #endif
 569        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 570        do i=1,nres
 571          write (iout,'(a3,i5,2f10.1)')
 572      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 573        enddo
 574 #ifdef MP
 575       endif
 576 #endif
 577       nct=nres
 578 cd      print *,'NNT=',NNT,' NCT=',NCT
 579       if (itype(1).eq.21) nnt=2
 580       if (itype(nres).eq.21) nct=nct-1
 581
 582 C     Juyong:READ init_vars
 583 C     Initialize variables!
 584 C     Juyong:READ read_info
 585 C     READ fragment information!!
 586 C     both routines should be in dfa.F file!!
 587
 588       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
 589      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
 590        call init_dfa_vars
 591        print*, 'init_dfa_vars finished!'
 592        call read_dfa_info
 593        print*, 'read_dfa_info finished!'
 594       endif
 595 C
 596 C
 597
 598
 599       if (pdbref) then
 600         if(me.eq.king.or..not.out1file)
 601      &   write (iout,'(a,i3)') 'nsup=',nsup
 602         nstart_seq=nnt
 603         if (nsup.le.(nct-nnt+1)) then
 604           do i=0,nct-nnt+1-nsup
 605             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 606               nstart_seq=nnt+i
 607               goto 111
 608             endif
 609           enddo
 610           write (iout,'(a)')
 611      &            'Error - sequences to be superposed do not match.'
 612           stop
 613         else
 614           do i=0,nsup-(nct-nnt+1)
 615             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 616      &      then
 617               nstart_sup=nstart_sup+i
 618               nsup=nct-nnt+1
 619               goto 111
 620             endif
 621           enddo
 622           write (iout,'(a)')
 623      &            'Error - sequences to be superposed do not match.'
 624         endif
 625   111   continue
 626         if (nsup.eq.0) nsup=nct-nnt
 627         if (nstart_sup.eq.0) nstart_sup=nnt
 628         if (nstart_seq.eq.0) nstart_seq=nnt
 629         if(me.eq.king.or..not.out1file)
 630      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 631      &                 ' nstart_seq=',nstart_seq
 632       endif
 633 c--- Zscore rms -------
 634       if (nz_start.eq.0) nz_start=nnt
 635       if (nz_end.eq.0 .and. nsup.gt.0) then
 636         nz_end=nnt+nsup-1
 637       else if (nz_end.eq.0) then
 638         nz_end=nct
 639       endif
 640       if(me.eq.king.or..not.out1file)then
 641        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 642        write (iout,*) 'IZ_SC=',iz_sc
 643       endif
 644 c----------------------
 645       call init_int_table
 646       if (refstr) then
 647         if (.not.pdbref) then
 648           call read_angles(inp,*38)
 649           goto 39
 650    38     write (iout,'(a)') 'Error reading reference structure.'
 651 #ifdef MPI
 652           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 653           stop 'Error reading reference structure'
 654 #endif
 655    39     call chainbuild
 656           call setup_var
 657 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 658           nstart_sup=nnt
 659           nstart_seq=nnt
 660           nsup=nct-nnt+1
 661           do i=1,2*nres
 662             do j=1,3
 663               cref(j,i)=c(j,i)
 664             enddo
 665           enddo
 666           call contact(.true.,ncont_ref,icont_ref,co)
 667         endif
 668 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 669         call flush(iout)
 670         if (constr_dist.gt.0) call read_dist_constr
 671 c        write (iout,*) "After read_dist_constr nhpb",nhpb
 672         call hpb_partition
 673         if(me.eq.king.or..not.out1file)
 674      &   write (iout,*) 'Contact order:',co
 675         if (pdbref) then
 676         if(me.eq.king.or..not.out1file)
 677      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 678         do i=1,ncont_ref
 679           do j=1,2
 680             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 681           enddo
 682           if(me.eq.king.or..not.out1file)
 683      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 684      &     icont_ref(1,i),' ',
 685      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 686         enddo
 687         endif
 688       endif
 689       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 690      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 691      &    modecalc.ne.10) then
 692 C If input structure hasn't been supplied from the PDB file read or generate
 693 C initial geometry.
 694         if (iranconf.eq.0 .and. .not. extconf) then
 695           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 696      &     write (iout,'(a)') 'Initial geometry will be read in.'
 697           if (read_cart) then
 698             read(inp,'(8f10.5)',end=36,err=36)
 699      &       ((c(l,k),l=1,3),k=1,nres),
 700      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 701             call int_from_cart1(.false.)
 702             do i=1,nres-1
 703               do j=1,3
 704                 dc(j,i)=c(j,i+1)-c(j,i)
 705                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 706               enddo
 707             enddo
 708             do i=nnt,nct
 709               if (itype(i).ne.10) then
 710                 do j=1,3
 711                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 712                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 713                 enddo
 714               endif
 715             enddo
 716             return
 717           else
 718             call read_angles(inp,*36)
 719           endif
 720           goto 37
 721    36     write (iout,'(a)') 'Error reading angle file.'
 722 #ifdef MPI
 723           call mpi_finalize( MPI_COMM_WORLD,IERR )
 724 #endif
 725           stop 'Error reading angle file.'
 726    37     continue
 727         else if (extconf) then
 728          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 729      &    write (iout,'(a)') 'Extended chain initial geometry.'
 730          do i=3,nres
 731           theta(i)=90d0*deg2rad
 732          enddo
 733          do i=4,nres
 734           phi(i)=180d0*deg2rad
 735          enddo
 736          do i=2,nres-1
 737           alph(i)=110d0*deg2rad
 738          enddo
 739          do i=2,nres-1
 740           omeg(i)=-120d0*deg2rad
 741           if (itype(i).le.0) omeg(i)=-omeg(i)
 742          enddo
 743         else
 744           if(me.eq.king.or..not.out1file)
 745      &     write (iout,'(a)') 'Random-generated initial geometry.'
 746
 747
 748 #ifdef MPI
 749           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 750      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 751 #endif
 752             do itrial=1,100
 753               itmp=1
 754 c              call gen_rand_conf(itmp,*30)
 755               goto 40
 756    30         write (iout,*) 'Failed to generate random conformation',
 757      &          ', itrial=',itrial
 758               write (*,*) 'Processor:',me,
 759      &          ' Failed to generate random conformation',
 760      &          ' itrial=',itrial
 761               call intout
 762
 763 #ifdef AIX
 764               call flush_(iout)
 765 #else
 766               call flush(iout)
 767 #endif
 768             enddo
 769             write (iout,'(a,i3,a)') 'Processor:',me,
 770      &        ' error in generating random conformation.'
 771             write (*,'(a,i3,a)') 'Processor:',me,
 772      &        ' error in generating random conformation.'
 773             call flush(iout)
 774 #ifdef MPI
 775             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 776    40       continue
 777           endif
 778 #else
 779           do itrial=1,100
 780             itmp=1
 781 c            call gen_rand_conf(itmp,*31)
 782             goto 40
 783    31       write (iout,*) 'Failed to generate random conformation',
 784      &        ', itrial=',itrial
 785             write (*,*) 'Failed to generate random conformation',
 786      &        ', itrial=',itrial
 787           enddo
 788           write (iout,'(a,i3,a)') 'Processor:',me,
 789      &      ' error in generating random conformation.'
 790           write (*,'(a,i3,a)') 'Processor:',me,
 791      &      ' error in generating random conformation.'
 792           stop
 793    40     continue
 794 #endif
 795         endif
 796       elseif (modecalc.eq.4) then
 797         read (inp,'(a)') intinname
 798         open (intin,file=intinname,status='old',err=333)
 799         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 800      &  write (iout,'(a)') 'intinname',intinname
 801         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 802         goto 334
 803   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 804 #ifdef MPI
 805         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 806 #endif
 807         stop 'Error opening angle file.'
 808   334   continue
 809
 810       endif
 811 C Generate distance constraints, if the PDB structure is to be regularized.
 812       if (nthread.gt.0) then
 813         call read_threadbase
 814       endif
 815       call setup_var
 816       if (me.eq.king .or. .not. out1file)
 817      & call intout
 818       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 819         write (iout,'(/a,i3,a)')
 820      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 821         write (iout,'(20i4)') (iss(i),i=1,ns)
 822         write (iout,'(/a/)') 'Pre-formed links are:'
 823         do i=1,nss
 824           i1=ihpb(i)-nres
 825           i2=jhpb(i)-nres
 826           it1=itype(i1)
 827           it2=itype(i2)
 828           if (me.eq.king.or..not.out1file)
 829      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 830      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 831      &    ebr,forcon(i)
 832         enddo
 833         write (iout,'(a)')
 834       endif
 835 c       if (i2ndstr.gt.0) call secstrp2dihc
 836 c      call geom_to_var(nvar,x)
 837 c      call etotal(energia(0))
 838 c      call enerprint(energia(0))
 839 c      call briefout(0,etot)
 840 c      stop
 841 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 842 cd    write (iout,'(a)') 'Variable list:'
 843 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 844 #ifdef MPI
 845       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 846      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 847      &  'Processor',myrank,': end reading molecular data.'
 848 #endif
 849       return
 850       end
 851 c--------------------------------------------------------------------------
 852       logical function seq_comp(itypea,itypeb,length)
 853       implicit none
 854       integer length,itypea(length),itypeb(length)
 855       integer i
 856       do i=1,length
 857         if (itypea(i).ne.itypeb(i)) then
 858           seq_comp=.false.
 859           return
 860         endif
 861       enddo
 862       seq_comp=.true.
 863       return
 864       end
 865 c-----------------------------------------------------------------------------
 866       subroutine read_bridge
 867 C Read information about disulfide bridges.
 868       implicit real*8 (a-h,o-z)
 869       include 'DIMENSIONS'
 870 #ifdef MPI
 871       include 'mpif.h'
 872 #endif
 873       include 'COMMON.IOUNITS'
 874       include 'COMMON.GEO'
 875       include 'COMMON.VAR'
 876       include 'COMMON.INTERACT'
 877       include 'COMMON.LOCAL'
 878       include 'COMMON.NAMES'
 879       include 'COMMON.CHAIN'
 880       include 'COMMON.FFIELD'
 881       include 'COMMON.SBRIDGE'
 882       include 'COMMON.HEADER'
 883       include 'COMMON.CONTROL'
 884 c      include 'COMMON.DBASE'
 885 c      include 'COMMON.THREAD'
 886       include 'COMMON.TIME1'
 887       include 'COMMON.SETUP'
 888 C Read bridging residues.
 889       read (inp,*) ns,(iss(i),i=1,ns)
 890       print *,'ns=',ns
 891       if(me.eq.king.or..not.out1file)
 892      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 893 C Check whether the specified bridging residues are cystines.
 894       do i=1,ns
 895         if (itype(iss(i)).ne.1) then
 896           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 897      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 898      &   ' can form a disulfide bridge?!!!'
 899           write (*,'(2a,i3,a)')
 900      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 901      &   ' can form a disulfide bridge?!!!'
 902 #ifdef MPI
 903          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 904          stop
 905 #endif
 906         endif
 907       enddo
 908 C Read preformed bridges.
 909       if (ns.gt.0) then
 910       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 911       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 912       if (nss.gt.0) then
 913         nhpb=nss
 914 C Check if the residues involved in bridges are in the specified list of
 915 C bridging residues.
 916         do i=1,nss
 917           do j=1,i-1
 918             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 919      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 920               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 921      &      ' contains residues present in other pairs.'
 922               write (*,'(a,i3,a)') 'Disulfide pair',i,
 923      &      ' contains residues present in other pairs.'
 924 #ifdef MPI
 925               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 926               stop
 927 #endif
 928             endif
 929           enddo
 930           do j=1,ns
 931             if (ihpb(i).eq.iss(j)) goto 10
 932           enddo
 933           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 934    10     continue
 935           do j=1,ns
 936             if (jhpb(i).eq.iss(j)) goto 20
 937           enddo
 938           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 939    20     continue
 940           dhpb(i)=dbr
 941           forcon(i)=fbr
 942         enddo
 943         do i=1,nss
 944           ihpb(i)=ihpb(i)+nres
 945           jhpb(i)=jhpb(i)+nres
 946         enddo
 947       endif
 948       endif
 949       return
 950       end
 951 c----------------------------------------------------------------------------
 952       subroutine read_x(kanal,*)
 953       implicit real*8 (a-h,o-z)
 954       include 'DIMENSIONS'
 955       include 'COMMON.GEO'
 956       include 'COMMON.VAR'
 957       include 'COMMON.CHAIN'
 958       include 'COMMON.IOUNITS'
 959       include 'COMMON.CONTROL'
 960       include 'COMMON.LOCAL'
 961       include 'COMMON.INTERACT'
 962 c Read coordinates from input
 963 c
 964       read(kanal,'(8f10.5)',end=10,err=10)
 965      &  ((c(l,k),l=1,3),k=1,nres),
 966      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 967       do j=1,3
 968         c(j,nres+1)=c(j,1)
 969         c(j,2*nres)=c(j,nres)
 970       enddo
 971       call int_from_cart1(.false.)
 972       do i=1,nres-1
 973         do j=1,3
 974           dc(j,i)=c(j,i+1)-c(j,i)
 975           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 976         enddo
 977       enddo
 978       do i=nnt,nct
 979         if (itype(i).ne.10) then
 980           do j=1,3
 981             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 982             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 983           enddo
 984         endif
 985       enddo
 986
 987       return
 988    10 return1
 989       end
 990 c------------------------------------------------------------------------------
 991       subroutine setup_var
 992       implicit real*8 (a-h,o-z)
 993       include 'DIMENSIONS'
 994       include 'COMMON.IOUNITS'
 995       include 'COMMON.GEO'
 996       include 'COMMON.VAR'
 997       include 'COMMON.INTERACT'
 998       include 'COMMON.LOCAL'
 999       include 'COMMON.NAMES'
1000       include 'COMMON.CHAIN'
1001       include 'COMMON.FFIELD'
1002       include 'COMMON.SBRIDGE'
1003       include 'COMMON.HEADER'
1004       include 'COMMON.CONTROL'
1005 c      include 'COMMON.DBASE'
1006 c      include 'COMMON.THREAD'
1007       include 'COMMON.TIME1'
1008 C Set up variable list.
1009       ntheta=nres-2
1010       nphi=nres-3
1011       nvar=ntheta+nphi
1012       nside=0
1013       do i=2,nres-1
1014         if (itype(i).ne.10) then
1015           nside=nside+1
1016           ialph(i,1)=nvar+nside
1017           ialph(nside,2)=i
1018         endif
1019       enddo
1020       if (indphi.gt.0) then
1021         nvar=nphi
1022       else if (indback.gt.0) then
1023         nvar=nphi+ntheta
1024       else
1025         nvar=nvar+2*nside
1026       endif
1027 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1028       return
1029       end
1030 c----------------------------------------------------------------------------
1031       subroutine gen_dist_constr
1032 C Generate CA distance constraints.
1033       implicit real*8 (a-h,o-z)
1034       include 'DIMENSIONS'
1035       include 'COMMON.IOUNITS'
1036       include 'COMMON.GEO'
1037       include 'COMMON.VAR'
1038       include 'COMMON.INTERACT'
1039       include 'COMMON.LOCAL'
1040       include 'COMMON.NAMES'
1041       include 'COMMON.CHAIN'
1042       include 'COMMON.FFIELD'
1043       include 'COMMON.SBRIDGE'
1044       include 'COMMON.HEADER'
1045       include 'COMMON.CONTROL'
1046 c      include 'COMMON.DBASE'
1047 c      include 'COMMON.THREAD'
1048       include 'COMMON.TIME1'
1049       dimension itype_pdb(maxres)
1050       common /pizda/ itype_pdb
1051       character*2 iden
1052 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1053 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1054 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1055 cd     & ' nsup',nsup
1056       do i=nstart_sup,nstart_sup+nsup-1
1057 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1058 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1059         do j=i+2,nstart_sup+nsup-1
1060           nhpb=nhpb+1
1061           ihpb(nhpb)=i+nstart_seq-nstart_sup
1062           jhpb(nhpb)=j+nstart_seq-nstart_sup
1063           forcon(nhpb)=weidis
1064           dhpb(nhpb)=dist(i,j)
1065         enddo
1066       enddo
1067 cd      write (iout,'(a)') 'Distance constraints:'
1068 cd      do i=nss+1,nhpb
1069 cd        ii=ihpb(i)
1070 cd        jj=jhpb(i)
1071 cd        iden='CA'
1072 cd        if (ii.gt.nres) then
1073 cd          iden='SC'
1074 cd          ii=ii-nres
1075 cd          jj=jj-nres
1076 cd        endif
1077 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1078 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1079 cd     &  dhpb(i),forcon(i)
1080 cd      enddo
1081       return
1082       end
1083 c----------------------------------------------------------------------------
1084       subroutine csaread
1085       implicit real*8 (a-h,o-z)
1086       include 'DIMENSIONS'
1087       include 'COMMON.IOUNITS'
1088       include 'COMMON.GEO'
1089       include 'COMMON.CSA'
1090       include 'COMMON.BANK'
1091       include 'COMMON.CONTROL'
1092       character*80 ucase
1093       character*620 mcmcard
1094       call card_concat(mcmcard)
1095
1096       call readi(mcmcard,'NCONF',nconf,50)
1097       call readi(mcmcard,'NADD',nadd,0)
1098       call readi(mcmcard,'JSTART',jstart,1)
1099       call readi(mcmcard,'JEND',jend,1)
1100       call readi(mcmcard,'NSTMAX',nstmax,500000)
1101       call readi(mcmcard,'N0',n0,1)
1102       call readi(mcmcard,'N1',n1,6)
1103       call readi(mcmcard,'N2',n2,4)
1104       call readi(mcmcard,'N3',n3,0)
1105       call readi(mcmcard,'N4',n4,0)
1106       call readi(mcmcard,'N5',n5,0)
1107       call readi(mcmcard,'N6',n6,10)
1108       call readi(mcmcard,'N7',n7,0)
1109       call readi(mcmcard,'N8',n8,0)
1110       call readi(mcmcard,'N9',n9,0)
1111       call readi(mcmcard,'N14',n14,0)
1112       call readi(mcmcard,'N15',n15,0)
1113       call readi(mcmcard,'N16',n16,0)
1114       call readi(mcmcard,'N17',n17,0)
1115       call readi(mcmcard,'N18',n18,0)
1116
1117       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1118
1119       call readi(mcmcard,'NDIFF',ndiff,2)
1120       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1121       call readi(mcmcard,'IS1',is1,1)
1122       call readi(mcmcard,'IS2',is2,8)
1123       call readi(mcmcard,'NRAN0',nran0,4)
1124       call readi(mcmcard,'NRAN1',nran1,2)
1125       call readi(mcmcard,'IRR',irr,1)
1126       call readi(mcmcard,'NSEED',nseed,20)
1127       call readi(mcmcard,'NTOTAL',ntotal,10000)
1128       call reada(mcmcard,'CUT1',cut1,2.0d0)
1129       call reada(mcmcard,'CUT2',cut2,5.0d0)
1130       call reada(mcmcard,'ESTOP',estop,-300000.0d0)
1131       call readi(mcmcard,'ICMAX',icmax,1)
1132       call readi(mcmcard,'NBANKM',nbankm,400)
1133       call readi(mcmcard,'IUCUT',iucut,2)
1134       call readi(mcmcard,'IRESTART',irestart,0)
1135 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1136       ntbankm=0
1137 c!bankt
1138       call reada(mcmcard,'DELE',dele,20.0d0)
1139       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1140       call readi(mcmcard,'IREF',iref,0)
1141       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1142       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1143       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1144       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1145       write (iout,*) "NCONF_IN",nconf_in
1146       tm_score=(index(mcmcard,'TMSCORE').gt.0)
1147       if (tm_score) write (iout,*) "Using TM_Score instead of DIFF",
1148      &      " for torsional angles"
1149       return
1150       end
1151
1152       subroutine read_minim
1153       implicit real*8 (a-h,o-z)
1154       include 'DIMENSIONS'
1155       include 'COMMON.MINIM'
1156       include 'COMMON.IOUNITS'
1157       character*80 ucase
1158       character*320 minimcard
1159       call card_concat(minimcard)
1160       call readi(minimcard,'MAXMIN',maxmin,2000)
1161       call readi(minimcard,'MAXFUN',maxfun,5000)
1162       call readi(minimcard,'MINMIN',minmin,maxmin)
1163       call readi(minimcard,'MINFUN',minfun,maxmin)
1164       call reada(minimcard,'TOLF',tolf,1.0D-2)
1165       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1166       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1167       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1168       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1169       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1170      &         'Options in energy minimization:'
1171       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1172      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1173      & 'MinMin:',MinMin,' MinFun:',MinFun,
1174      & ' TolF:',TolF,' RTolF:',RTolF
1175       return
1176       end
1177 c----------------------------------------------------------------------------
1178       subroutine read_angles(kanal,*)
1179       implicit real*8 (a-h,o-z)
1180       include 'DIMENSIONS'
1181       include 'COMMON.GEO'
1182       include 'COMMON.VAR'
1183       include 'COMMON.CHAIN'
1184       include 'COMMON.IOUNITS'
1185       include 'COMMON.CONTROL'
1186 c Read angles from input
1187 c
1188        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1189        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1190        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1191        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1192
1193        do i=1,nres
1194 c 9/7/01 avoid 180 deg valence angle
1195         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1196 c
1197         theta(i)=deg2rad*theta(i)
1198         phi(i)=deg2rad*phi(i)
1199         alph(i)=deg2rad*alph(i)
1200         omeg(i)=deg2rad*omeg(i)
1201        enddo
1202       return
1203    10 return1
1204       end
1205 c----------------------------------------------------------------------------
1206       subroutine reada(rekord,lancuch,wartosc,default)
1207       implicit none
1208       character*(*) rekord,lancuch
1209       double precision wartosc,default
1210       integer ilen,iread
1211       external ilen
1212       iread=index(rekord,lancuch)
1213       if (iread.eq.0) then
1214         wartosc=default
1215         return
1216       endif
1217       iread=iread+ilen(lancuch)+1
1218       read (rekord(iread:),*,err=10,end=10) wartosc
1219       return
1220   10  wartosc=default
1221       return
1222       end
1223 c----------------------------------------------------------------------------
1224       subroutine readi(rekord,lancuch,wartosc,default)
1225       implicit none
1226       character*(*) rekord,lancuch
1227       integer wartosc,default
1228       integer ilen,iread
1229       external ilen
1230       iread=index(rekord,lancuch)
1231       if (iread.eq.0) then
1232         wartosc=default
1233         return
1234       endif
1235       iread=iread+ilen(lancuch)+1
1236       read (rekord(iread:),*,err=10,end=10) wartosc
1237       return
1238   10  wartosc=default
1239       return
1240       end
1241 c----------------------------------------------------------------------------
1242       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1243       implicit none
1244       integer dim,i
1245       integer tablica(dim),default
1246       character*(*) rekord,lancuch
1247       character*80 aux
1248       integer ilen,iread
1249       external ilen
1250       do i=1,dim
1251         tablica(i)=default
1252       enddo
1253       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1254       if (iread.eq.0) return
1255       iread=iread+ilen(lancuch)+1
1256       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1257    10 return
1258       end
1259 c----------------------------------------------------------------------------
1260       subroutine multreada(rekord,lancuch,tablica,dim,default)
1261       implicit none
1262       integer dim,i
1263       double precision tablica(dim),default
1264       character*(*) rekord,lancuch
1265       character*80 aux
1266       integer ilen,iread
1267       external ilen
1268       do i=1,dim
1269         tablica(i)=default
1270       enddo
1271       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1272       if (iread.eq.0) return
1273       iread=iread+ilen(lancuch)+1
1274       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1275    10 return
1276       end
1277 c----------------------------------------------------------------------------
1278       subroutine openunits
1279       implicit real*8 (a-h,o-z)
1280       include 'DIMENSIONS'
1281 #ifdef MPI
1282       include 'mpif.h'
1283       character*16 form,nodename
1284       integer nodelen
1285 #endif
1286       include 'COMMON.SETUP'
1287       include 'COMMON.IOUNITS'
1288 c      include 'COMMON.MD'
1289       include 'COMMON.CONTROL'
1290       integer lenpre,lenpot,ilen,lentmp
1291       external ilen
1292       character*3 out1file_text,ucase
1293       character*3 ll
1294       external ucase
1295 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1296       call getenv_loc("PREFIX",prefix)
1297       pref_orig = prefix
1298       call getenv_loc("POT",pot)
1299       call getenv_loc("DIRTMP",tmpdir)
1300       call getenv_loc("CURDIR",curdir)
1301       call getenv_loc("OUT1FILE",out1file_text)
1302 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1303       out1file_text=ucase(out1file_text)
1304       if (out1file_text(1:1).eq."Y") then
1305         out1file=.true.
1306       else
1307         out1file=fg_rank.gt.0
1308       endif
1309       lenpre=ilen(prefix)
1310       lenpot=ilen(pot)
1311       lentmp=ilen(tmpdir)
1312       if (lentmp.gt.0) then
1313           write (*,'(80(1h!))')
1314           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1315           write (*,'(80(1h!))')
1316           write (*,*)"All output files will be on node /tmp directory."
1317 #ifdef MPI
1318         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1319         if (me.eq.king) then
1320           write (*,*) "The master node is ",nodename
1321         else if (fg_rank.eq.0) then
1322           write (*,*) "I am the CG slave node ",nodename
1323         else
1324           write (*,*) "I am the FG slave node ",nodename
1325         endif
1326 #endif
1327         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1328         lenpre = lentmp+lenpre+1
1329       endif
1330       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1331 C Get the names and open the input files
1332 #if defined(WINIFL) || defined(WINPGI)
1333       open(1,file=pref_orig(:ilen(pref_orig))//
1334      &  '.inp',status='old',readonly,shared)
1335        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1336 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1337 C Get parameter filenames and open the parameter files.
1338       call getenv_loc('BONDPAR',bondname)
1339       open (ibond,file=bondname,status='old',readonly,shared)
1340       call getenv_loc('THETPAR',thetname)
1341       open (ithep,file=thetname,status='old',readonly,shared)
1342 #ifndef CRYST_THETA
1343       call getenv_loc('THETPARPDB',thetname_pdb)
1344       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1345 #endif
1346       call getenv_loc('ROTPAR',rotname)
1347       open (irotam,file=rotname,status='old',readonly,shared)
1348 #ifndef CRYST_SC
1349       call getenv_loc('ROTPARPDB',rotname_pdb)
1350       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1351 #endif
1352       call getenv_loc('TORPAR',torname)
1353       open (itorp,file=torname,status='old',readonly,shared)
1354       call getenv_loc('TORDPAR',tordname)
1355       open (itordp,file=tordname,status='old',readonly,shared)
1356       call getenv_loc('FOURIER',fouriername)
1357       open (ifourier,file=fouriername,status='old',readonly,shared)
1358       call getenv_loc('ELEPAR',elename)
1359       open (ielep,file=elename,status='old',readonly,shared)
1360       call getenv_loc('SIDEPAR',sidename)
1361       open (isidep,file=sidename,status='old',readonly,shared)
1362 #elif (defined CRAY) || (defined AIX)
1363       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1364      &  action='read')
1365 c      print *,"Processor",myrank," opened file 1"
1366       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1367 c      print *,"Processor",myrank," opened file 9"
1368 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1369 C Get parameter filenames and open the parameter files.
1370       call getenv_loc('BONDPAR',bondname)
1371       open (ibond,file=bondname,status='old',action='read')
1372 c      print *,"Processor",myrank," opened file IBOND"
1373       call getenv_loc('THETPAR',thetname)
1374       open (ithep,file=thetname,status='old',action='read')
1375 c      print *,"Processor",myrank," opened file ITHEP"
1376 #ifndef CRYST_THETA
1377       call getenv_loc('THETPARPDB',thetname_pdb)
1378       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1379 #endif
1380       call getenv_loc('ROTPAR',rotname)
1381       open (irotam,file=rotname,status='old',action='read')
1382 c      print *,"Processor",myrank," opened file IROTAM"
1383 #ifndef CRYST_SC
1384       call getenv_loc('ROTPARPDB',rotname_pdb)
1385       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1386 #endif
1387       call getenv_loc('TORPAR',torname)
1388       open (itorp,file=torname,status='old',action='read')
1389 c      print *,"Processor",myrank," opened file ITORP"
1390       call getenv_loc('TORDPAR',tordname)
1391       open (itordp,file=tordname,status='old',action='read')
1392 c      print *,"Processor",myrank," opened file ITORDP"
1393       call getenv_loc('SCCORPAR',sccorname)
1394       open (isccor,file=sccorname,status='old',action='read')
1395 c      print *,"Processor",myrank," opened file ISCCOR"
1396       call getenv_loc('FOURIER',fouriername)
1397       open (ifourier,file=fouriername,status='old',action='read')
1398 c      print *,"Processor",myrank," opened file IFOURIER"
1399       call getenv_loc('ELEPAR',elename)
1400       open (ielep,file=elename,status='old',action='read')
1401 c      print *,"Processor",myrank," opened file IELEP"
1402       call getenv_loc('SIDEPAR',sidename)
1403       open (isidep,file=sidename,status='old',action='read')
1404 c      print *,"Processor",myrank," opened file ISIDEP"
1405 c      print *,"Processor",myrank," opened parameter files"
1406 #elif (defined G77)
1407       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1408       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1409 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1410 C Get parameter filenames and open the parameter files.
1411       call getenv_loc('BONDPAR',bondname)
1412       open (ibond,file=bondname,status='old')
1413       call getenv_loc('THETPAR',thetname)
1414       open (ithep,file=thetname,status='old')
1415 #ifndef CRYST_THETA
1416       call getenv_loc('THETPARPDB',thetname_pdb)
1417       open (ithep_pdb,file=thetname_pdb,status='old')
1418 #endif
1419       call getenv_loc('ROTPAR',rotname)
1420       open (irotam,file=rotname,status='old')
1421 #ifndef CRYST_SC
1422       call getenv_loc('ROTPARPDB',rotname_pdb)
1423       open (irotam_pdb,file=rotname_pdb,status='old')
1424 #endif
1425       call getenv_loc('TORPAR',torname)
1426       open (itorp,file=torname,status='old')
1427       call getenv_loc('TORDPAR',tordname)
1428       open (itordp,file=tordname,status='old')
1429       call getenv_loc('SCCORPAR',sccorname)
1430       open (isccor,file=sccorname,status='old')
1431       call getenv_loc('FOURIER',fouriername)
1432       open (ifourier,file=fouriername,status='old')
1433       call getenv_loc('ELEPAR',elename)
1434       open (ielep,file=elename,status='old')
1435       call getenv_loc('SIDEPAR',sidename)
1436       open (isidep,file=sidename,status='old')
1437 #else
1438       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1439      &  readonly)
1440        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1441 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1442 C Get parameter filenames and open the parameter files.
1443       call getenv_loc('BONDPAR',bondname)
1444       open (ibond,file=bondname,status='old',readonly)
1445       call getenv_loc('THETPAR',thetname)
1446       open (ithep,file=thetname,status='old',readonly)
1447 #ifndef CRYST_THETA
1448       call getenv_loc('THETPARPDB',thetname_pdb)
1449       print *,"thetname_pdb ",thetname_pdb
1450       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1451       print *,ithep_pdb," opened"
1452 #endif
1453       call getenv_loc('ROTPAR',rotname)
1454       open (irotam,file=rotname,status='old',readonly)
1455 #ifndef CRYST_SC
1456       call getenv_loc('ROTPARPDB',rotname_pdb)
1457       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1458 #endif
1459       call getenv_loc('TORPAR',torname)
1460       open (itorp,file=torname,status='old',readonly)
1461       call getenv_loc('TORDPAR',tordname)
1462       open (itordp,file=tordname,status='old',readonly)
1463       call getenv_loc('SCCORPAR',sccorname)
1464       open (isccor,file=sccorname,status='old',readonly)
1465       call getenv_loc('FOURIER',fouriername)
1466       open (ifourier,file=fouriername,status='old',readonly)
1467       call getenv_loc('ELEPAR',elename)
1468       open (ielep,file=elename,status='old',readonly)
1469       call getenv_loc('SIDEPAR',sidename)
1470       open (isidep,file=sidename,status='old',readonly)
1471 #endif
1472 #ifndef OLDSCP
1473 C
1474 C 8/9/01 In the newest version SCp interaction constants are read from a file
1475 C Use -DOLDSCP to use hard-coded constants instead.
1476 C
1477       call getenv_loc('SCPPAR',scpname)
1478 #if defined(WINIFL) || defined(WINPGI)
1479       open (iscpp,file=scpname,status='old',readonly,shared)
1480 #elif (defined CRAY)  || (defined AIX)
1481       open (iscpp,file=scpname,status='old',action='read')
1482 #elif (defined G77)
1483       open (iscpp,file=scpname,status='old')
1484 #else
1485       open (iscpp,file=scpname,status='old',readonly)
1486 #endif
1487 #endif
1488       call getenv_loc('PATTERN',patname)
1489 #if defined(WINIFL) || defined(WINPGI)
1490       open (icbase,file=patname,status='old',readonly,shared)
1491 #elif (defined CRAY)  || (defined AIX)
1492       open (icbase,file=patname,status='old',action='read')
1493 #elif (defined G77)
1494       open (icbase,file=patname,status='old')
1495 #else
1496       open (icbase,file=patname,status='old',readonly)
1497 #endif
1498 #ifdef MPI
1499 C Open output file only for CG processes
1500 c      print *,"Processor",myrank," fg_rank",fg_rank
1501       if (fg_rank.eq.0) then
1502
1503       if (nodes.eq.1) then
1504         npos=3
1505       else
1506         npos = dlog10(dfloat(nodes-1))+1
1507       endif
1508       if (npos.lt.3) npos=3
1509       write (liczba,'(i1)') npos
1510       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1511      &  //')'
1512       write (liczba,form) me
1513       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1514      &  liczba(:ilen(liczba))
1515       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1516      &  //'.int'
1517       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1518      &  //'.pdb'
1519       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1520      &  liczba(:ilen(liczba))//'.mol2'
1521       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1522      &  liczba(:ilen(liczba))//'.stat'
1523       if (lentmp.gt.0)
1524      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1525      &      //liczba(:ilen(liczba))//'.stat')
1526       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1527      &  //'.rst'
1528 c      if(usampl) then
1529 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1530 c     & liczba(:ilen(liczba))//'.const'
1531 c      endif
1532
1533       endif
1534 #else
1535       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1536       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1537       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1538       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1539       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1540       if (lentmp.gt.0)
1541      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1542      &    '.stat')
1543       rest2name=prefix(:ilen(prefix))//'.rst'
1544 c      if(usampl) then
1545 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1546 c      endif
1547 #endif
1548 #if defined(AIX) || defined(PGI)
1549       if (me.eq.king .or. .not. out1file)
1550      &   open(iout,file=outname,status='unknown')
1551 #ifdef DEBUG
1552       if (fg_rank.gt.0) then
1553         write (liczba,'(i3.3)') myrank/nfgtasks
1554         write (ll,'(bz,i3.3)') fg_rank
1555         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1556      &   status='unknown')
1557       endif
1558 #endif
1559       if(me.eq.king) then
1560        open(igeom,file=intname,status='unknown',position='append')
1561        open(ipdb,file=pdbname,status='unknown')
1562        open(imol2,file=mol2name,status='unknown')
1563        open(istat,file=statname,status='unknown',position='append')
1564       else
1565 c1out       open(iout,file=outname,status='unknown')
1566       endif
1567 #else
1568       if (me.eq.king .or. .not.out1file)
1569      &    open(iout,file=outname,status='unknown')
1570 #ifdef DEBUG
1571       if (fg_rank.gt.0) then
1572         write (liczba,'(i3.3)') myrank/nfgtasks
1573         write (ll,'(bz,i3.3)') fg_rank
1574         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1575      &   status='unknown')
1576       endif
1577 #endif
1578       if(me.eq.king) then
1579        open(igeom,file=intname,status='unknown',access='append')
1580        open(ipdb,file=pdbname,status='unknown')
1581        open(imol2,file=mol2name,status='unknown')
1582        open(istat,file=statname,status='unknown',access='append')
1583       else
1584 c1out       open(iout,file=outname,status='unknown')
1585       endif
1586 #endif
1587       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1588       csa_seed=prefix(:lenpre)//'.CSA.seed'
1589       csa_history=prefix(:lenpre)//'.CSA.history'
1590       csa_bank=prefix(:lenpre)//'.CSA.bank'
1591       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1592       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1593       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1594 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1595       csa_int=prefix(:lenpre)//'.int'
1596       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1597       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1598       csa_in=prefix(:lenpre)//'.CSA.in'
1599 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1600 C Write file names
1601       if (me.eq.king)then
1602       write (iout,'(80(1h-))')
1603       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1604       write (iout,'(80(1h-))')
1605       write (iout,*) "Input file                      : ",
1606      &  pref_orig(:ilen(pref_orig))//'.inp'
1607       write (iout,*) "Output file                     : ",
1608      &  outname(:ilen(outname))
1609       write (iout,*)
1610       write (iout,*) "Sidechain potential file        : ",
1611      &  sidename(:ilen(sidename))
1612 #ifndef OLDSCP
1613       write (iout,*) "SCp potential file              : ",
1614      &  scpname(:ilen(scpname))
1615 #endif
1616       write (iout,*) "Electrostatic potential file    : ",
1617      &  elename(:ilen(elename))
1618       write (iout,*) "Cumulant coefficient file       : ",
1619      &  fouriername(:ilen(fouriername))
1620       write (iout,*) "Torsional parameter file        : ",
1621      &  torname(:ilen(torname))
1622       write (iout,*) "Double torsional parameter file : ",
1623      &  tordname(:ilen(tordname))
1624       write (iout,*) "SCCOR parameter file : ",
1625      &  sccorname(:ilen(sccorname))
1626       write (iout,*) "Bond & inertia constant file    : ",
1627      &  bondname(:ilen(bondname))
1628       write (iout,*) "Bending parameter file          : ",
1629      &  thetname(:ilen(thetname))
1630       write (iout,*) "Rotamer parameter file          : ",
1631      &  rotname(:ilen(rotname))
1632       write (iout,*) "Threading database              : ",
1633      &  patname(:ilen(patname))
1634       if (lentmp.ne.0)
1635      &write (iout,*)" DIRTMP                          : ",
1636      &  tmpdir(:lentmp)
1637       write (iout,'(80(1h-))')
1638       endif
1639       return
1640       end
1641 c----------------------------------------------------------------------------
1642       subroutine card_concat(card)
1643       implicit real*8 (a-h,o-z)
1644       include 'DIMENSIONS'
1645       include 'COMMON.IOUNITS'
1646       character*(*) card
1647       character*80 karta,ucase
1648       external ilen
1649       read (inp,'(a)') karta
1650       karta=ucase(karta)
1651       card=' '
1652       do while (karta(80:80).eq.'&')
1653         card=card(:ilen(card)+1)//karta(:79)
1654         read (inp,'(a)') karta
1655         karta=ucase(karta)
1656       enddo
1657       card=card(:ilen(card)+1)//karta
1658       return
1659       end
1660
1661       subroutine read_dist_constr
1662       implicit real*8 (a-h,o-z)
1663       include 'DIMENSIONS'
1664 #ifdef MPI
1665       include 'mpif.h'
1666 #endif
1667       include 'COMMON.SETUP'
1668       include 'COMMON.CONTROL'
1669       include 'COMMON.CHAIN'
1670       include 'COMMON.IOUNITS'
1671       include 'COMMON.SBRIDGE'
1672       integer ifrag_(2,100),ipair_(2,100)
1673       double precision wfrag_(100),wpair_(100)
1674       character*500 controlcard
1675       write (iout,*) "Calling read_dist_constr"
1676       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1677       call flush(iout)
1678       call card_concat(controlcard)
1679       call readi(controlcard,"NFRAG",nfrag_,0)
1680       call readi(controlcard,"NPAIR",npair_,0)
1681       call readi(controlcard,"NDIST",ndist_,0)
1682       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1683       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1684       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1685       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1686       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1687 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1688 c      write (iout,*) "IFRAG"
1689 c      do i=1,nfrag_
1690 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1691 c      enddo
1692 c      write (iout,*) "IPAIR"
1693 c      do i=1,npair_
1694 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1695 c      enddo
1696       call flush(iout)
1697       do i=1,nfrag_
1698         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1699         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1700      &    ifrag_(2,i)=nstart_sup+nsup-1
1701 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1702         call flush(iout)
1703         if (wfrag_(i).gt.0.0d0) then
1704         do j=ifrag_(1,i),ifrag_(2,i)-1
1705           do k=j+1,ifrag_(2,i)
1706             write (iout,*) "j",j," k",k
1707             ddjk=dist(j,k)
1708             if (constr_dist.eq.1) then
1709             nhpb=nhpb+1
1710             ihpb(nhpb)=j
1711             jhpb(nhpb)=k
1712               dhpb(nhpb)=ddjk
1713             forcon(nhpb)=wfrag_(i)
1714             else if (constr_dist.eq.2) then
1715               if (ddjk.le.dist_cut) then
1716                 nhpb=nhpb+1
1717                 ihpb(nhpb)=j
1718                 jhpb(nhpb)=k
1719                 dhpb(nhpb)=ddjk
1720                 forcon(nhpb)=wfrag_(i)
1721               endif
1722             else
1723               nhpb=nhpb+1
1724               ihpb(nhpb)=j
1725               jhpb(nhpb)=k
1726               dhpb(nhpb)=ddjk
1727               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1728             endif
1729 #ifdef MPI
1730             if (.not.out1file .or. me.eq.king)
1731      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1732      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1733 #else
1734             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1735      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1736 #endif
1737           enddo
1738         enddo
1739         endif
1740       enddo
1741       do i=1,npair_
1742         if (wpair_(i).gt.0.0d0) then
1743         ii = ipair_(1,i)
1744         jj = ipair_(2,i)
1745         if (ii.gt.jj) then
1746           itemp=ii
1747           ii=jj
1748           jj=itemp
1749         endif
1750         do j=ifrag_(1,ii),ifrag_(2,ii)
1751           do k=ifrag_(1,jj),ifrag_(2,jj)
1752             nhpb=nhpb+1
1753             ihpb(nhpb)=j
1754             jhpb(nhpb)=k
1755             forcon(nhpb)=wpair_(i)
1756             dhpb(nhpb)=dist(j,k)
1757 #ifdef MPI
1758             if (.not.out1file .or. me.eq.king)
1759      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1760      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1761 #else
1762             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1763      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1764 #endif
1765           enddo
1766         enddo
1767         endif
1768       enddo
1769       do i=1,ndist_
1770         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1771         if (forcon(nhpb+1).gt.0.0d0) then
1772           nhpb=nhpb+1
1773           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1774 #ifdef MPI
1775           if (.not.out1file .or. me.eq.king)
1776      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1777      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1778 #else
1779           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1780      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1781 #endif
1782         endif
1783       enddo
1784       call flush(iout)
1785       return
1786       end
1787 c-------------------------------------------------------------------------------
1788 #ifdef WINIFL
1789       subroutine flush(iu)
1790       return
1791       end
1792 #endif
1793 #ifdef AIX
1794       subroutine flush(iu)
1795       call flush_(iu)
1796       return
1797       end
1798 #endif
1799 c------------------------------------------------------------------------------
1800       subroutine copy_to_tmp(source)
1801       include "DIMENSIONS"
1802       include "COMMON.IOUNITS"
1803       character*(*) source
1804       character* 256 tmpfile
1805       integer ilen
1806       external ilen
1807       logical ex
1808       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1809       inquire(file=tmpfile,exist=ex)
1810       if (ex) then
1811         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1812      &   " to temporary directory..."
1813         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1814         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1815       endif
1816       return
1817       end
1818 c------------------------------------------------------------------------------
1819       subroutine move_from_tmp(source)
1820       include "DIMENSIONS"
1821       include "COMMON.IOUNITS"
1822       character*(*) source
1823       integer ilen
1824       external ilen
1825       write (*,*) "Moving ",source(:ilen(source)),
1826      & " from temporary directory to working directory"
1827       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1828       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1829       return
1830       end
1831 c------------------------------------------------------------------------------
1832       subroutine random_init(seed)
1833 C
1834 C Initialize random number generator
1835 C
1836       implicit real*8 (a-h,o-z)
1837       include 'DIMENSIONS'
1838 #ifdef AMD64
1839       integer*8 iseedi8
1840 #endif
1841 #ifdef MPI
1842       include 'mpif.h'
1843       logical OKRandom, prng_restart
1844       real*8  r1
1845       integer iseed_array(4)
1846 #endif
1847       include 'COMMON.IOUNITS'
1848       include 'COMMON.TIME1'
1849 c      include 'COMMON.THREAD'
1850       include 'COMMON.SBRIDGE'
1851       include 'COMMON.CONTROL'
1852       include 'COMMON.MCM'
1853 c      include 'COMMON.MAP'
1854       include 'COMMON.HEADER'
1855 c      include 'COMMON.CSA'
1856       include 'COMMON.CHAIN'
1857 c      include 'COMMON.MUCA'
1858 c      include 'COMMON.MD'
1859       include 'COMMON.FFIELD'
1860       include 'COMMON.SETUP'
1861       iseed=-dint(dabs(seed))
1862       if (iseed.eq.0) then
1863         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1864      &    'Random seed undefined. The program will stop.'
1865         write (*,'(/80(1h*)/20x,a/80(1h*))')
1866      &    'Random seed undefined. The program will stop.'
1867 #ifdef MPI
1868         call mpi_finalize(mpi_comm_world,ierr)
1869 #endif
1870         stop 'Bad random seed.'
1871       endif
1872 #ifdef MPI
1873       if (fg_rank.eq.0) then
1874       seed=seed*(me+1)+1
1875 #ifdef AMD64
1876       iseedi8=dint(seed)
1877       if(me.eq.king .or. .not. out1file)
1878      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1879       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1880       OKRandom = prng_restart(me,iseedi8)
1881 #else
1882       do i=1,4
1883        tmp=65536.0d0**(4-i)
1884        iseed_array(i) = dint(seed/tmp)
1885        seed=seed-iseed_array(i)*tmp
1886       enddo
1887       if(me.eq.king .or. .not. out1file)
1888      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1889      &                 (iseed_array(i),i=1,4)
1890       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1891      &                 (iseed_array(i),i=1,4)
1892       OKRandom = prng_restart(me,iseed_array)
1893 #endif
1894       if (OKRandom) then
1895         r1=ran_number(0.0D0,1.0D0)
1896         if(me.eq.king .or. .not. out1file)
1897      &   write (iout,*) 'ran_num',r1
1898         if (r1.lt.0.0d0) OKRandom=.false.
1899       endif
1900       if (.not.OKRandom) then
1901         write (iout,*) 'PRNG IS NOT WORKING!!!'
1902         print *,'PRNG IS NOT WORKING!!!'
1903         if (me.eq.0) then
1904          call flush(iout)
1905          call mpi_abort(mpi_comm_world,error_msg,ierr)
1906          stop
1907         else
1908          write (iout,*) 'too many processors for parallel prng'
1909          write (*,*) 'too many processors for parallel prng'
1910          call flush(iout)
1911          stop
1912         endif
1913       endif
1914       endif
1915 #else
1916       call vrndst(iseed)
1917       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1918 #endif
1919       return
1920       end