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