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       return
1135       end
1136
1137       subroutine read_minim
1138       implicit real*8 (a-h,o-z)
1139       include 'DIMENSIONS'
1140       include 'COMMON.MINIM'
1141       include 'COMMON.IOUNITS'
1142       character*80 ucase
1143       character*320 minimcard
1144       call card_concat(minimcard)
1145       call readi(minimcard,'MAXMIN',maxmin,2000)
1146       call readi(minimcard,'MAXFUN',maxfun,5000)
1147       call readi(minimcard,'MINMIN',minmin,maxmin)
1148       call readi(minimcard,'MINFUN',minfun,maxmin)
1149       call reada(minimcard,'TOLF',tolf,1.0D-2)
1150       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1151       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1152      &         'Options in energy minimization:'
1153       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1154      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1155      & 'MinMin:',MinMin,' MinFun:',MinFun,
1156      & ' TolF:',TolF,' RTolF:',RTolF
1157       return
1158       end
1159 c----------------------------------------------------------------------------
1160       subroutine read_angles(kanal,*)
1161       implicit real*8 (a-h,o-z)
1162       include 'DIMENSIONS'
1163       include 'COMMON.GEO'
1164       include 'COMMON.VAR'
1165       include 'COMMON.CHAIN'
1166       include 'COMMON.IOUNITS'
1167       include 'COMMON.CONTROL'
1168 c Read angles from input
1169 c
1170        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1171        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1172        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1173        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1174
1175        do i=1,nres
1176 c 9/7/01 avoid 180 deg valence angle
1177         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1178 c
1179         theta(i)=deg2rad*theta(i)
1180         phi(i)=deg2rad*phi(i)
1181         alph(i)=deg2rad*alph(i)
1182         omeg(i)=deg2rad*omeg(i)
1183        enddo
1184       return
1185    10 return1
1186       end
1187 c----------------------------------------------------------------------------
1188       subroutine reada(rekord,lancuch,wartosc,default)
1189       implicit none
1190       character*(*) rekord,lancuch
1191       double precision wartosc,default
1192       integer ilen,iread
1193       external ilen
1194       iread=index(rekord,lancuch)
1195       if (iread.eq.0) then
1196         wartosc=default
1197         return
1198       endif
1199       iread=iread+ilen(lancuch)+1
1200       read (rekord(iread:),*,err=10,end=10) wartosc
1201       return
1202   10  wartosc=default
1203       return
1204       end
1205 c----------------------------------------------------------------------------
1206       subroutine readi(rekord,lancuch,wartosc,default)
1207       implicit none
1208       character*(*) rekord,lancuch
1209       integer wartosc,default
1210       integer ilen,iread
1211       external ilen
1212       iread=index(rekord,lancuch)
1213       if (iread.eq.0) then
1214         wartosc=default
1215         return
1216       endif
1217       iread=iread+ilen(lancuch)+1
1218       read (rekord(iread:),*,err=10,end=10) wartosc
1219       return
1220   10  wartosc=default
1221       return
1222       end
1223 c----------------------------------------------------------------------------
1224       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1225       implicit none
1226       integer dim,i
1227       integer tablica(dim),default
1228       character*(*) rekord,lancuch
1229       character*80 aux
1230       integer ilen,iread
1231       external ilen
1232       do i=1,dim
1233         tablica(i)=default
1234       enddo
1235       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1236       if (iread.eq.0) return
1237       iread=iread+ilen(lancuch)+1
1238       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1239    10 return
1240       end
1241 c----------------------------------------------------------------------------
1242       subroutine multreada(rekord,lancuch,tablica,dim,default)
1243       implicit none
1244       integer dim,i
1245       double precision tablica(dim),default
1246       character*(*) rekord,lancuch
1247       character*80 aux
1248       integer ilen,iread
1249       external ilen
1250       do i=1,dim
1251         tablica(i)=default
1252       enddo
1253       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1254       if (iread.eq.0) return
1255       iread=iread+ilen(lancuch)+1
1256       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1257    10 return
1258       end
1259 c----------------------------------------------------------------------------
1260       subroutine openunits
1261       implicit real*8 (a-h,o-z)
1262       include 'DIMENSIONS'
1263 #ifdef MPI
1264       include 'mpif.h'
1265       character*16 form,nodename
1266       integer nodelen
1267 #endif
1268       include 'COMMON.SETUP'
1269       include 'COMMON.IOUNITS'
1270 c      include 'COMMON.MD'
1271       include 'COMMON.CONTROL'
1272       integer lenpre,lenpot,ilen,lentmp
1273       external ilen
1274       character*3 out1file_text,ucase
1275       character*3 ll
1276       external ucase
1277 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1278       call getenv_loc("PREFIX",prefix)
1279       pref_orig = prefix
1280       call getenv_loc("POT",pot)
1281       call getenv_loc("DIRTMP",tmpdir)
1282       call getenv_loc("CURDIR",curdir)
1283       call getenv_loc("OUT1FILE",out1file_text)
1284 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1285       out1file_text=ucase(out1file_text)
1286       if (out1file_text(1:1).eq."Y") then
1287         out1file=.true.
1288       else
1289         out1file=fg_rank.gt.0
1290       endif
1291       lenpre=ilen(prefix)
1292       lenpot=ilen(pot)
1293       lentmp=ilen(tmpdir)
1294       if (lentmp.gt.0) then
1295           write (*,'(80(1h!))')
1296           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1297           write (*,'(80(1h!))')
1298           write (*,*)"All output files will be on node /tmp directory."
1299 #ifdef MPI
1300         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1301         if (me.eq.king) then
1302           write (*,*) "The master node is ",nodename
1303         else if (fg_rank.eq.0) then
1304           write (*,*) "I am the CG slave node ",nodename
1305         else
1306           write (*,*) "I am the FG slave node ",nodename
1307         endif
1308 #endif
1309         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1310         lenpre = lentmp+lenpre+1
1311       endif
1312       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1313 C Get the names and open the input files
1314 #if defined(WINIFL) || defined(WINPGI)
1315       open(1,file=pref_orig(:ilen(pref_orig))//
1316      &  '.inp',status='old',readonly,shared)
1317        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1318 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1319 C Get parameter filenames and open the parameter files.
1320       call getenv_loc('BONDPAR',bondname)
1321       open (ibond,file=bondname,status='old',readonly,shared)
1322       call getenv_loc('THETPAR',thetname)
1323       open (ithep,file=thetname,status='old',readonly,shared)
1324 #ifndef CRYST_THETA
1325       call getenv_loc('THETPARPDB',thetname_pdb)
1326       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1327 #endif
1328       call getenv_loc('ROTPAR',rotname)
1329       open (irotam,file=rotname,status='old',readonly,shared)
1330 #ifndef CRYST_SC
1331       call getenv_loc('ROTPARPDB',rotname_pdb)
1332       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
1333 #endif
1334       call getenv_loc('TORPAR',torname)
1335       open (itorp,file=torname,status='old',readonly,shared)
1336       call getenv_loc('TORDPAR',tordname)
1337       open (itordp,file=tordname,status='old',readonly,shared)
1338       call getenv_loc('FOURIER',fouriername)
1339       open (ifourier,file=fouriername,status='old',readonly,shared)
1340       call getenv_loc('ELEPAR',elename)
1341       open (ielep,file=elename,status='old',readonly,shared)
1342       call getenv_loc('SIDEPAR',sidename)
1343       open (isidep,file=sidename,status='old',readonly,shared)
1344 #elif (defined CRAY) || (defined AIX)
1345       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1346      &  action='read')
1347 c      print *,"Processor",myrank," opened file 1"
1348       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1349 c      print *,"Processor",myrank," opened file 9"
1350 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1351 C Get parameter filenames and open the parameter files.
1352       call getenv_loc('BONDPAR',bondname)
1353       open (ibond,file=bondname,status='old',action='read')
1354 c      print *,"Processor",myrank," opened file IBOND"
1355       call getenv_loc('THETPAR',thetname)
1356       open (ithep,file=thetname,status='old',action='read')
1357 c      print *,"Processor",myrank," opened file ITHEP"
1358 #ifndef CRYST_THETA
1359       call getenv_loc('THETPARPDB',thetname_pdb)
1360       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1361 #endif
1362       call getenv_loc('ROTPAR',rotname)
1363       open (irotam,file=rotname,status='old',action='read')
1364 c      print *,"Processor",myrank," opened file IROTAM"
1365 #ifndef CRYST_SC
1366       call getenv_loc('ROTPARPDB',rotname_pdb)
1367       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1368 #endif
1369       call getenv_loc('TORPAR',torname)
1370       open (itorp,file=torname,status='old',action='read')
1371 c      print *,"Processor",myrank," opened file ITORP"
1372       call getenv_loc('TORDPAR',tordname)
1373       open (itordp,file=tordname,status='old',action='read')
1374 c      print *,"Processor",myrank," opened file ITORDP"
1375       call getenv_loc('SCCORPAR',sccorname)
1376       open (isccor,file=sccorname,status='old',action='read')
1377 c      print *,"Processor",myrank," opened file ISCCOR"
1378       call getenv_loc('FOURIER',fouriername)
1379       open (ifourier,file=fouriername,status='old',action='read')
1380 c      print *,"Processor",myrank," opened file IFOURIER"
1381       call getenv_loc('ELEPAR',elename)
1382       open (ielep,file=elename,status='old',action='read')
1383 c      print *,"Processor",myrank," opened file IELEP"
1384       call getenv_loc('SIDEPAR',sidename)
1385       open (isidep,file=sidename,status='old',action='read')
1386 c      print *,"Processor",myrank," opened file ISIDEP"
1387 c      print *,"Processor",myrank," opened parameter files"
1388 #elif (defined G77)
1389       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1390       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1391 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1392 C Get parameter filenames and open the parameter files.
1393       call getenv_loc('BONDPAR',bondname)
1394       open (ibond,file=bondname,status='old')
1395       call getenv_loc('THETPAR',thetname)
1396       open (ithep,file=thetname,status='old')
1397 #ifndef CRYST_THETA
1398       call getenv_loc('THETPARPDB',thetname_pdb)
1399       open (ithep_pdb,file=thetname_pdb,status='old')
1400 #endif
1401       call getenv_loc('ROTPAR',rotname)
1402       open (irotam,file=rotname,status='old')
1403 #ifndef CRYST_SC
1404       call getenv_loc('ROTPARPDB',rotname_pdb)
1405       open (irotam_pdb,file=rotname_pdb,status='old')
1406 #endif
1407       call getenv_loc('TORPAR',torname)
1408       open (itorp,file=torname,status='old')
1409       call getenv_loc('TORDPAR',tordname)
1410       open (itordp,file=tordname,status='old')
1411       call getenv_loc('SCCORPAR',sccorname)
1412       open (isccor,file=sccorname,status='old')
1413       call getenv_loc('FOURIER',fouriername)
1414       open (ifourier,file=fouriername,status='old')
1415       call getenv_loc('ELEPAR',elename)
1416       open (ielep,file=elename,status='old')
1417       call getenv_loc('SIDEPAR',sidename)
1418       open (isidep,file=sidename,status='old')
1419 #else
1420       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1421      &  readonly)
1422        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1423 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1424 C Get parameter filenames and open the parameter files.
1425       call getenv_loc('BONDPAR',bondname)
1426       open (ibond,file=bondname,status='old',readonly)
1427       call getenv_loc('THETPAR',thetname)
1428       open (ithep,file=thetname,status='old',readonly)
1429 #ifndef CRYST_THETA
1430       call getenv_loc('THETPARPDB',thetname_pdb)
1431       print *,"thetname_pdb ",thetname_pdb
1432       open (ithep_pdb,file=thetname_pdb,status='old',readonly)
1433       print *,ithep_pdb," opened"
1434 #endif
1435       call getenv_loc('ROTPAR',rotname)
1436       open (irotam,file=rotname,status='old',readonly)
1437 #ifndef CRYST_SC
1438       call getenv_loc('ROTPARPDB',rotname_pdb)
1439       open (irotam_pdb,file=rotname_pdb,status='old',readonly)
1440 #endif
1441       call getenv_loc('TORPAR',torname)
1442       open (itorp,file=torname,status='old',readonly)
1443       call getenv_loc('TORDPAR',tordname)
1444       open (itordp,file=tordname,status='old',readonly)
1445       call getenv_loc('SCCORPAR',sccorname)
1446       open (isccor,file=sccorname,status='old',readonly)
1447       call getenv_loc('FOURIER',fouriername)
1448       open (ifourier,file=fouriername,status='old',readonly)
1449       call getenv_loc('ELEPAR',elename)
1450       open (ielep,file=elename,status='old',readonly)
1451       call getenv_loc('SIDEPAR',sidename)
1452       open (isidep,file=sidename,status='old',readonly)
1453 #endif
1454 #ifndef OLDSCP
1455 C
1456 C 8/9/01 In the newest version SCp interaction constants are read from a file
1457 C Use -DOLDSCP to use hard-coded constants instead.
1458 C
1459       call getenv_loc('SCPPAR',scpname)
1460 #if defined(WINIFL) || defined(WINPGI)
1461       open (iscpp,file=scpname,status='old',readonly,shared)
1462 #elif (defined CRAY)  || (defined AIX)
1463       open (iscpp,file=scpname,status='old',action='read')
1464 #elif (defined G77)
1465       open (iscpp,file=scpname,status='old')
1466 #else
1467       open (iscpp,file=scpname,status='old',readonly)
1468 #endif
1469 #endif
1470       call getenv_loc('PATTERN',patname)
1471 #if defined(WINIFL) || defined(WINPGI)
1472       open (icbase,file=patname,status='old',readonly,shared)
1473 #elif (defined CRAY)  || (defined AIX)
1474       open (icbase,file=patname,status='old',action='read')
1475 #elif (defined G77)
1476       open (icbase,file=patname,status='old')
1477 #else
1478       open (icbase,file=patname,status='old',readonly)
1479 #endif
1480 #ifdef MPI
1481 C Open output file only for CG processes
1482 c      print *,"Processor",myrank," fg_rank",fg_rank
1483       if (fg_rank.eq.0) then
1484
1485       if (nodes.eq.1) then
1486         npos=3
1487       else
1488         npos = dlog10(dfloat(nodes-1))+1
1489       endif
1490       if (npos.lt.3) npos=3
1491       write (liczba,'(i1)') npos
1492       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
1493      &  //')'
1494       write (liczba,form) me
1495       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
1496      &  liczba(:ilen(liczba))
1497       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1498      &  //'.int'
1499       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
1500      &  //'.pdb'
1501       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
1502      &  liczba(:ilen(liczba))//'.mol2'
1503       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
1504      &  liczba(:ilen(liczba))//'.stat'
1505       if (lentmp.gt.0)
1506      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
1507      &      //liczba(:ilen(liczba))//'.stat')
1508       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
1509      &  //'.rst'
1510 c      if(usampl) then
1511 c          qname=prefix(:lenpre)//'_'//pot(:lenpot)//
1512 c     & liczba(:ilen(liczba))//'.const'
1513 c      endif
1514
1515       endif
1516 #else
1517       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
1518       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
1519       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
1520       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
1521       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
1522       if (lentmp.gt.0)
1523      & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)//
1524      &    '.stat')
1525       rest2name=prefix(:ilen(prefix))//'.rst'
1526 c      if(usampl) then
1527 c         qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
1528 c      endif
1529 #endif
1530 #if defined(AIX) || defined(PGI)
1531       if (me.eq.king .or. .not. out1file)
1532      &   open(iout,file=outname,status='unknown')
1533 #ifdef DEBUG
1534       if (fg_rank.gt.0) then
1535         write (liczba,'(i3.3)') myrank/nfgtasks
1536         write (ll,'(bz,i3.3)') fg_rank
1537         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1538      &   status='unknown')
1539       endif
1540 #endif
1541       if(me.eq.king) then
1542        open(igeom,file=intname,status='unknown',position='append')
1543        open(ipdb,file=pdbname,status='unknown')
1544        open(imol2,file=mol2name,status='unknown')
1545        open(istat,file=statname,status='unknown',position='append')
1546       else
1547 c1out       open(iout,file=outname,status='unknown')
1548       endif
1549 #else
1550       if (me.eq.king .or. .not.out1file)
1551      &    open(iout,file=outname,status='unknown')
1552 #ifdef DEBUG
1553       if (fg_rank.gt.0) then
1554         write (liczba,'(i3.3)') myrank/nfgtasks
1555         write (ll,'(bz,i3.3)') fg_rank
1556         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
1557      &   status='unknown')
1558       endif
1559 #endif
1560       if(me.eq.king) then
1561        open(igeom,file=intname,status='unknown',access='append')
1562        open(ipdb,file=pdbname,status='unknown')
1563        open(imol2,file=mol2name,status='unknown')
1564        open(istat,file=statname,status='unknown',access='append')
1565       else
1566 c1out       open(iout,file=outname,status='unknown')
1567       endif
1568 #endif
1569       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
1570       csa_seed=prefix(:lenpre)//'.CSA.seed'
1571       csa_history=prefix(:lenpre)//'.CSA.history'
1572       csa_bank=prefix(:lenpre)//'.CSA.bank'
1573       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
1574       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
1575       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
1576 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
1577       csa_int=prefix(:lenpre)//'.int'
1578       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
1579       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
1580       csa_in=prefix(:lenpre)//'.CSA.in'
1581 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
1582 C Write file names
1583       if (me.eq.king)then
1584       write (iout,'(80(1h-))')
1585       write (iout,'(30x,a)') "FILE ASSIGNMENT"
1586       write (iout,'(80(1h-))')
1587       write (iout,*) "Input file                      : ",
1588      &  pref_orig(:ilen(pref_orig))//'.inp'
1589       write (iout,*) "Output file                     : ",
1590      &  outname(:ilen(outname))
1591       write (iout,*)
1592       write (iout,*) "Sidechain potential file        : ",
1593      &  sidename(:ilen(sidename))
1594 #ifndef OLDSCP
1595       write (iout,*) "SCp potential file              : ",
1596      &  scpname(:ilen(scpname))
1597 #endif
1598       write (iout,*) "Electrostatic potential file    : ",
1599      &  elename(:ilen(elename))
1600       write (iout,*) "Cumulant coefficient file       : ",
1601      &  fouriername(:ilen(fouriername))
1602       write (iout,*) "Torsional parameter file        : ",
1603      &  torname(:ilen(torname))
1604       write (iout,*) "Double torsional parameter file : ",
1605      &  tordname(:ilen(tordname))
1606       write (iout,*) "SCCOR parameter file : ",
1607      &  sccorname(:ilen(sccorname))
1608       write (iout,*) "Bond & inertia constant file    : ",
1609      &  bondname(:ilen(bondname))
1610       write (iout,*) "Bending parameter file          : ",
1611      &  thetname(:ilen(thetname))
1612       write (iout,*) "Rotamer parameter file          : ",
1613      &  rotname(:ilen(rotname))
1614       write (iout,*) "Threading database              : ",
1615      &  patname(:ilen(patname))
1616       if (lentmp.ne.0)
1617      &write (iout,*)" DIRTMP                          : ",
1618      &  tmpdir(:lentmp)
1619       write (iout,'(80(1h-))')
1620       endif
1621       return
1622       end
1623 c----------------------------------------------------------------------------
1624       subroutine card_concat(card)
1625       implicit real*8 (a-h,o-z)
1626       include 'DIMENSIONS'
1627       include 'COMMON.IOUNITS'
1628       character*(*) card
1629       character*80 karta,ucase
1630       external ilen
1631       read (inp,'(a)') karta
1632       karta=ucase(karta)
1633       card=' '
1634       do while (karta(80:80).eq.'&')
1635         card=card(:ilen(card)+1)//karta(:79)
1636         read (inp,'(a)') karta
1637         karta=ucase(karta)
1638       enddo
1639       card=card(:ilen(card)+1)//karta
1640       return
1641       end
1642
1643       subroutine read_dist_constr
1644       implicit real*8 (a-h,o-z)
1645       include 'DIMENSIONS'
1646 #ifdef MPI
1647       include 'mpif.h'
1648 #endif
1649       include 'COMMON.SETUP'
1650       include 'COMMON.CONTROL'
1651       include 'COMMON.CHAIN'
1652       include 'COMMON.IOUNITS'
1653       include 'COMMON.SBRIDGE'
1654       integer ifrag_(2,100),ipair_(2,100)
1655       double precision wfrag_(100),wpair_(100)
1656       character*500 controlcard
1657       write (iout,*) "Calling read_dist_constr"
1658       write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1659       call flush(iout)
1660       call card_concat(controlcard)
1661       call readi(controlcard,"NFRAG",nfrag_,0)
1662       call readi(controlcard,"NPAIR",npair_,0)
1663       call readi(controlcard,"NDIST",ndist_,0)
1664       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1665       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1666       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1667       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1668       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1669 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1670 c      write (iout,*) "IFRAG"
1671 c      do i=1,nfrag_
1672 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1673 c      enddo
1674 c      write (iout,*) "IPAIR"
1675 c      do i=1,npair_
1676 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1677 c      enddo
1678       call flush(iout)
1679       do i=1,nfrag_
1680         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1681         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1682      &    ifrag_(2,i)=nstart_sup+nsup-1
1683 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1684         call flush(iout)
1685         if (wfrag_(i).gt.0.0d0) then
1686         do j=ifrag_(1,i),ifrag_(2,i)-1
1687           do k=j+1,ifrag_(2,i)
1688             write (iout,*) "j",j," k",k
1689             ddjk=dist(j,k)
1690             if (constr_dist.eq.1) then
1691             nhpb=nhpb+1
1692             ihpb(nhpb)=j
1693             jhpb(nhpb)=k
1694               dhpb(nhpb)=ddjk
1695             forcon(nhpb)=wfrag_(i)
1696             else if (constr_dist.eq.2) then
1697               if (ddjk.le.dist_cut) then
1698                 nhpb=nhpb+1
1699                 ihpb(nhpb)=j
1700                 jhpb(nhpb)=k
1701                 dhpb(nhpb)=ddjk
1702                 forcon(nhpb)=wfrag_(i)
1703               endif
1704             else
1705               nhpb=nhpb+1
1706               ihpb(nhpb)=j
1707               jhpb(nhpb)=k
1708               dhpb(nhpb)=ddjk
1709               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1710             endif
1711 #ifdef MPI
1712             if (.not.out1file .or. me.eq.king)
1713      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1714      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1715 #else
1716             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1717      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1718 #endif
1719           enddo
1720         enddo
1721         endif
1722       enddo
1723       do i=1,npair_
1724         if (wpair_(i).gt.0.0d0) then
1725         ii = ipair_(1,i)
1726         jj = ipair_(2,i)
1727         if (ii.gt.jj) then
1728           itemp=ii
1729           ii=jj
1730           jj=itemp
1731         endif
1732         do j=ifrag_(1,ii),ifrag_(2,ii)
1733           do k=ifrag_(1,jj),ifrag_(2,jj)
1734             nhpb=nhpb+1
1735             ihpb(nhpb)=j
1736             jhpb(nhpb)=k
1737             forcon(nhpb)=wpair_(i)
1738             dhpb(nhpb)=dist(j,k)
1739 #ifdef MPI
1740             if (.not.out1file .or. me.eq.king)
1741      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1742      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1743 #else
1744             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1745      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1746 #endif
1747           enddo
1748         enddo
1749         endif
1750       enddo
1751       do i=1,ndist_
1752         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
1753         if (forcon(nhpb+1).gt.0.0d0) then
1754           nhpb=nhpb+1
1755           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1756 #ifdef MPI
1757           if (.not.out1file .or. me.eq.king)
1758      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1759      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1760 #else
1761           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1762      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1763 #endif
1764         endif
1765       enddo
1766       call flush(iout)
1767       return
1768       end
1769 c-------------------------------------------------------------------------------
1770 #ifdef WINIFL
1771       subroutine flush(iu)
1772       return
1773       end
1774 #endif
1775 #ifdef AIX
1776       subroutine flush(iu)
1777       call flush_(iu)
1778       return
1779       end
1780 #endif
1781 c------------------------------------------------------------------------------
1782       subroutine copy_to_tmp(source)
1783       include "DIMENSIONS"
1784       include "COMMON.IOUNITS"
1785       character*(*) source
1786       character* 256 tmpfile
1787       integer ilen
1788       external ilen
1789       logical ex
1790       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1791       inquire(file=tmpfile,exist=ex)
1792       if (ex) then
1793         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1794      &   " to temporary directory..."
1795         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1796         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1797       endif
1798       return
1799       end
1800 c------------------------------------------------------------------------------
1801       subroutine move_from_tmp(source)
1802       include "DIMENSIONS"
1803       include "COMMON.IOUNITS"
1804       character*(*) source
1805       integer ilen
1806       external ilen
1807       write (*,*) "Moving ",source(:ilen(source)),
1808      & " from temporary directory to working directory"
1809       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1810       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1811       return
1812       end
1813 c------------------------------------------------------------------------------
1814       subroutine random_init(seed)
1815 C
1816 C Initialize random number generator
1817 C
1818       implicit real*8 (a-h,o-z)
1819       include 'DIMENSIONS'
1820 #ifdef AMD64
1821       integer*8 iseedi8
1822 #endif
1823 #ifdef MPI
1824       include 'mpif.h'
1825       logical OKRandom, prng_restart
1826       real*8  r1
1827       integer iseed_array(4)
1828 #endif
1829       include 'COMMON.IOUNITS'
1830       include 'COMMON.TIME1'
1831 c      include 'COMMON.THREAD'
1832       include 'COMMON.SBRIDGE'
1833       include 'COMMON.CONTROL'
1834       include 'COMMON.MCM'
1835 c      include 'COMMON.MAP'
1836       include 'COMMON.HEADER'
1837 c      include 'COMMON.CSA'
1838       include 'COMMON.CHAIN'
1839 c      include 'COMMON.MUCA'
1840 c      include 'COMMON.MD'
1841       include 'COMMON.FFIELD'
1842       include 'COMMON.SETUP'
1843       iseed=-dint(dabs(seed))
1844       if (iseed.eq.0) then
1845         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1846      &    'Random seed undefined. The program will stop.'
1847         write (*,'(/80(1h*)/20x,a/80(1h*))')
1848      &    'Random seed undefined. The program will stop.'
1849 #ifdef MPI
1850         call mpi_finalize(mpi_comm_world,ierr)
1851 #endif
1852         stop 'Bad random seed.'
1853       endif
1854 #ifdef MPI
1855       if (fg_rank.eq.0) then
1856       seed=seed*(me+1)+1
1857 #ifdef AMD64
1858       iseedi8=dint(seed)
1859       if(me.eq.king .or. .not. out1file)
1860      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1861       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1862       OKRandom = prng_restart(me,iseedi8)
1863 #else
1864       do i=1,4
1865        tmp=65536.0d0**(4-i)
1866        iseed_array(i) = dint(seed/tmp)
1867        seed=seed-iseed_array(i)*tmp
1868       enddo
1869       if(me.eq.king .or. .not. out1file)
1870      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1871      &                 (iseed_array(i),i=1,4)
1872       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1873      &                 (iseed_array(i),i=1,4)
1874       OKRandom = prng_restart(me,iseed_array)
1875 #endif
1876       if (OKRandom) then
1877         r1=ran_number(0.0D0,1.0D0)
1878         if(me.eq.king .or. .not. out1file)
1879      &   write (iout,*) 'ran_num',r1
1880         if (r1.lt.0.0d0) OKRandom=.false.
1881       endif
1882       if (.not.OKRandom) then
1883         write (iout,*) 'PRNG IS NOT WORKING!!!'
1884         print *,'PRNG IS NOT WORKING!!!'
1885         if (me.eq.0) then
1886          call flush(iout)
1887          call mpi_abort(mpi_comm_world,error_msg,ierr)
1888          stop
1889         else
1890          write (iout,*) 'too many processors for parallel prng'
1891          write (*,*) 'too many processors for parallel prng'
1892          call flush(iout)
1893          stop
1894         endif
1895       endif
1896       endif
1897 #else
1898       call vrndst(iseed)
1899       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1900 #endif
1901       return
1902       end