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