source/unres/src_Eshel/readrtns.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 molecule information, molecule geometry, energy-term weights, and
  19 C restraints if requested
  20       call molread
  21 C Print restraint information
  22 #ifdef MPI
  23       if (.not. out1file .or. me.eq.king) then
  24 #endif
  25       if (nhpb.gt.nss)
  26      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  27      & " distance constraints have been imposed"
  28       do i=nss+1,nhpb
  29         write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
  30      &     ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
  31       enddo
  32 #ifdef MPI
  33       endif
  34 #endif
  35 c      print *,"Processor",myrank," leaves READRTNS"
  36       return
  37       end
  38 C-------------------------------------------------------------------------------
  39       subroutine read_control
  40 C
  41 C Read contorl data
  42 C
  43       implicit real*8 (a-h,o-z)
  44       include 'DIMENSIONS'
  45 #ifdef MP
  46       include 'mpif.h'
  47       logical OKRandom, prng_restart
  48       real*8  r1
  49 #endif
  50       include 'COMMON.IOUNITS'
  51       include 'COMMON.TIME1'
  52       include 'COMMON.THREAD'
  53       include 'COMMON.SBRIDGE'
  54       include 'COMMON.CONTROL'
  55       include 'COMMON.MCM'
  56       include 'COMMON.MAP'
  57       include 'COMMON.HEADER'
  58 csa      include 'COMMON.CSA'
  59       include 'COMMON.CHAIN'
  60       include 'COMMON.MUCA'
  61       include 'COMMON.MD'
  62       include 'COMMON.FFIELD'
  63       include 'COMMON.SETUP'
  64       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  65       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  66       character*80 ucase
  67       character*320 controlcard
  68
  69       nglob_csa=0
  70       eglob_csa=1d99
  71       nmin_csa=0
  72       read (INP,'(a)') titel
  73       call card_concat(controlcard)
  74 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  75 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  76       call reada(controlcard,'SEED',seed,0.0D0)
  77       call random_init(seed)
  78 C Set up the time limit (caution! The time must be input in minutes!)
  79       read_cart=index(controlcard,'READ_CART').gt.0
  80       catrace=index(controlcard,'CATRACE').gt.0
  81       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  82       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
  83       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
  84       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
  85       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
  86       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
  87       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
  88       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
  89       call reada(controlcard,'DRMS',drms,0.1D0)
  90       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
  91        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
  92        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
  93        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
  94        write (iout,'(a,f10.1)')'DRMS    = ',drms
  95        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
  96        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
  97       endif
  98       call readi(controlcard,'NZ_START',nz_start,0)
  99       call readi(controlcard,'NZ_END',nz_end,0)
 100       call readi(controlcard,'IZ_SC',iz_sc,0)
 101       timlim=60.0D0*timlim
 102       safety = 60.0d0*safety
 103       timem=timlim
 104       modecalc=0
 105       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 106       minim=(index(controlcard,'MINIMIZE').gt.0)
 107       dccart=(index(controlcard,'CART').gt.0)
 108       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 109       overlapsc=.not.overlapsc
 110       searchsc=(index(controlcard,'SEARCHSC').gt.0)
 111       sideadd=(index(controlcard,'SIDEADD').gt.0) .or. catrace
 112       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 113       outpdb=(index(controlcard,'PDBOUT').gt.0)
 114       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 115       pdbref=(index(controlcard,'PDBREF').gt.0)
 116       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 117       indpdb=index(controlcard,'PDBSTART')
 118       extconf=(index(controlcard,'EXTCONF').gt.0)
 119       call readi(controlcard,'IPRINT',iprint,0)
 120       call readi(controlcard,'MAXGEN',maxgen,10000)
 121       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 122       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 123       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 124      & write (iout,*) "RESCALE_MODE",rescale_mode
 125       if (index(controlcard,'REGULAR').gt.0.0D0) then
 126         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 127         modecalc=1
 128         refstr=.true.
 129         regular=.true.
 130       endif
 131       if (index(controlcard,'CHECKGRAD').gt.0) then
 132         modecalc=5
 133         if (index(controlcard,'CART').gt.0) then
 134           icheckgrad=1
 135         elseif (index(controlcard,'CARINT').gt.0) then
 136           icheckgrad=2
 137         else
 138           icheckgrad=3
 139         endif
 140       elseif (index(controlcard,'THREAD').gt.0) then
 141         modecalc=2
 142         call readi(controlcard,'THREAD',nthread,0)
 143         if (nthread.gt.0) then
 144           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 145         else
 146           if (fg_rank.eq.0)
 147      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 148           stop 'Error termination in Read_Control.'
 149         endif
 150       else if (index(controlcard,'MCMA').gt.0) then
 151         modecalc=3
 152       else if (index(controlcard,'MCEE').gt.0) then
 153         modecalc=6
 154       else if (index(controlcard,'MULTCONF').gt.0) then
 155         modecalc=4
 156       else if (index(controlcard,'MAP').gt.0) then
 157         modecalc=7
 158         call readi(controlcard,'MAP',nmap,0)
 159       else if (index(controlcard,'CSA').gt.0) then
 160            write(*,*) "CSA not supported in this version"
 161            stop
 162       else if (index(controlcard,'SOFTREG').gt.0) then
 163         modecalc=11
 164       else if (index(controlcard,'CHECK_BOND').gt.0) then
 165         modecalc=-1
 166       else if (index(controlcard,'TEST').gt.0) then
 167         modecalc=-2
 168       endif
 169
 170       lmuca=index(controlcard,'MUCA').gt.0
 171       call readi(controlcard,'MUCADYN',mucadyn,0)
 172       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 173       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 174      & then
 175        write (iout,*) 'MUCADYN=',mucadyn
 176        write (iout,*) 'MUCASMOOTH=',muca_smooth
 177       endif
 178
 179       iscode=index(controlcard,'ONE_LETTER')
 180       indphi=index(controlcard,'PHI')
 181       indback=index(controlcard,'BACK')
 182       iranconf=index(controlcard,'RAND_CONF')
 183       i2ndstr=index(controlcard,'USE_SEC_PRED')
 184       gradout=index(controlcard,'GRADOUT').gt.0
 185       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 186
 187       if(me.eq.king.or..not.out1file)
 188      & write (iout,'(2a)') diagmeth(kdiag),
 189      &  ' routine used to diagonalize matrices.'
 190       return
 191       end
 192 c--------------------------------------------------------------------------
 193       subroutine molread
 194 C
 195 C Read molecular data.
 196 C
 197       implicit real*8 (a-h,o-z)
 198       include 'DIMENSIONS'
 199 #ifdef MPI
 200       include 'mpif.h'
 201       integer error_msg
 202 #endif
 203       include 'COMMON.IOUNITS'
 204       include 'COMMON.GEO'
 205       include 'COMMON.VAR'
 206       include 'COMMON.INTERACT'
 207       include 'COMMON.LOCAL'
 208       include 'COMMON.NAMES'
 209       include 'COMMON.CHAIN'
 210       include 'COMMON.FFIELD'
 211       include 'COMMON.SBRIDGE'
 212       include 'COMMON.HEADER'
 213       include 'COMMON.CONTROL'
 214       include 'COMMON.DBASE'
 215       include 'COMMON.THREAD'
 216       include 'COMMON.CONTACTS'
 217       include 'COMMON.TORCNSTR'
 218       include 'COMMON.TIME1'
 219       include 'COMMON.BOUNDS'
 220       include 'COMMON.MD'
 221       include 'COMMON.REMD'
 222       include 'COMMON.SETUP'
 223       character*4 sequence(maxres)
 224       integer rescode
 225       double precision x(maxvar)
 226       character*320 weightcard
 227       character*80 weightcard_t,ucase
 228       dimension itype_pdb(maxres)
 229       common /pizda/ itype_pdb
 230       logical seq_comp,fail
 231       double precision energia(0:n_ene)
 232       integer ilen
 233       external ilen
 234 C
 235 C Body
 236 C
 237 C Read weights of the subsequent energy terms.
 238       if(hremd.gt.0) then
 239
 240        k=0
 241        do il=1,hremd
 242         do i=1,nrep
 243          do j=1,remd_m(i)
 244           i2set(k)=il
 245           k=k+1
 246          enddo
 247         enddo
 248        enddo
 249
 250        if(me.eq.king.or..not.out1file) then
 251         write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
 252         write (iout,*) 'Current weights for processor ',
 253      &                 me,' set ',i2set(me)
 254        endif
 255
 256        do i=1,hremd
 257          call card_concat(weightcard)
 258          call reada(weightcard,'WLONG',wlong,1.0D0)
 259          call reada(weightcard,'WSC',wsc,wlong)
 260          call reada(weightcard,'WSCP',wscp,wlong)
 261          call reada(weightcard,'WELEC',welec,1.0D0)
 262          call reada(weightcard,'WVDWPP',wvdwpp,welec)
 263          call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 264          call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 265          call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 266          call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 267          call reada(weightcard,'WTURN3',wturn3,1.0D0)
 268          call reada(weightcard,'WTURN4',wturn4,1.0D0)
 269          call reada(weightcard,'WTURN6',wturn6,1.0D0)
 270          call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 271          call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 272          call reada(weightcard,'WBOND',wbond,1.0D0)
 273          call reada(weightcard,'WTOR',wtor,1.0D0)
 274          call reada(weightcard,'WTORD',wtor_d,1.0D0)
 275          call reada(weightcard,'WANG',wang,1.0D0)
 276          call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 277          call reada(weightcard,'SCAL14',scal14,0.4D0)
 278          call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 279          call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 280          call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 281          call reada(weightcard,'TEMP0',temp0,300.0d0)
 282          if (index(weightcard,'SOFT').gt.0) ipot=6
 283 C 12/1/95 Added weight for the multi-body term WCORR
 284          call reada(weightcard,'WCORRH',wcorr,1.0D0)
 285          if (wcorr4.gt.0.0d0) wcorr=wcorr4
 286
 287          hweights(i,1)=wsc
 288          hweights(i,2)=wscp
 289          hweights(i,3)=welec
 290          hweights(i,4)=wcorr
 291          hweights(i,5)=wcorr5
 292          hweights(i,6)=wcorr6
 293          hweights(i,7)=wel_loc
 294          hweights(i,8)=wturn3
 295          hweights(i,9)=wturn4
 296          hweights(i,10)=wturn6
 297          hweights(i,11)=wang
 298          hweights(i,12)=wscloc
 299          hweights(i,13)=wtor
 300          hweights(i,14)=wtor_d
 301          hweights(i,15)=wstrain
 302          hweights(i,16)=wvdwpp
 303          hweights(i,17)=wbond
 304          hweights(i,18)=scal14
 305          hweights(i,21)=wsccor
 306
 307        enddo
 308
 309        do i=1,n_ene
 310          weights(i)=hweights(i2set(me),i)
 311        enddo
 312        wsc    =weights(1)
 313        wscp   =weights(2)
 314        welec  =weights(3)
 315        wcorr  =weights(4)
 316        wcorr5 =weights(5)
 317        wcorr6 =weights(6)
 318        wel_loc=weights(7)
 319        wturn3 =weights(8)
 320        wturn4 =weights(9)
 321        wturn6 =weights(10)
 322        wang   =weights(11)
 323        wscloc =weights(12)
 324        wtor   =weights(13)
 325        wtor_d =weights(14)
 326        wstrain=weights(15)
 327        wvdwpp =weights(16)
 328        wbond  =weights(17)
 329        scal14 =weights(18)
 330        wsccor =weights(21)
 331
 332
 333       else
 334        call card_concat(weightcard)
 335        call reada(weightcard,'WLONG',wlong,1.0D0)
 336        call reada(weightcard,'WSC',wsc,wlong)
 337        call reada(weightcard,'WSCP',wscp,wlong)
 338        call reada(weightcard,'WELEC',welec,1.0D0)
 339        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 340        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 341        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 342        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 343        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 344        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 345        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 346        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 347        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 348        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 349        call reada(weightcard,'WBOND',wbond,1.0D0)
 350        call reada(weightcard,'WTOR',wtor,1.0D0)
 351        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 352        call reada(weightcard,'WANG',wang,1.0D0)
 353        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 354        call reada(weightcard,'SCAL14',scal14,0.4D0)
 355        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 356        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 357        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 358        call reada(weightcard,'TEMP0',temp0,300.0d0)
 359        if (index(weightcard,'SOFT').gt.0) ipot=6
 360 C 12/1/95 Added weight for the multi-body term WCORR
 361        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 362        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 363        weights(1)=wsc
 364        weights(2)=wscp
 365        weights(3)=welec
 366        weights(4)=wcorr
 367        weights(5)=wcorr5
 368        weights(6)=wcorr6
 369        weights(7)=wel_loc
 370        weights(8)=wturn3
 371        weights(9)=wturn4
 372        weights(10)=wturn6
 373        weights(11)=wang
 374        weights(12)=wscloc
 375        weights(13)=wtor
 376        weights(14)=wtor_d
 377        weights(15)=wstrain
 378        weights(16)=wvdwpp
 379        weights(17)=wbond
 380        weights(18)=scal14
 381        weights(21)=wsccor
 382       endif
 383
 384       if(me.eq.king.or..not.out1file)
 385      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 386      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 387      &  wturn4,wturn6
 388    10 format (/'Energy-term weights (unscaled):'//
 389      & 'WSCC=   ',f10.6,' (SC-SC)'/
 390      & 'WSCP=   ',f10.6,' (SC-p)'/
 391      & 'WELEC=  ',f10.6,' (p-p electr)'/
 392      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 393      & 'WBOND=  ',f10.6,' (stretching)'/
 394      & 'WANG=   ',f10.6,' (bending)'/
 395      & 'WSCLOC= ',f10.6,' (SC local)'/
 396      & 'WTOR=   ',f10.6,' (torsional)'/
 397      & 'WTORD=  ',f10.6,' (double torsional)'/
 398      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 399      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 400      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 401      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 402      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 403      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 404      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 405      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 406      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 407       if(me.eq.king.or..not.out1file)then
 408        if (wcorr4.gt.0.0d0) then
 409         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 410      &   'between contact pairs of peptide groups'
 411         write (iout,'(2(a,f5.3/))')
 412      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 413      &  'Range of quenching the correlation terms:',2*delt_corr
 414        else if (wcorr.gt.0.0d0) then
 415         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 416      &   'between contact pairs of peptide groups'
 417        endif
 418        write (iout,'(a,f8.3)')
 419      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 420        write (iout,'(a,f8.3)')
 421      &  'General scaling factor of SC-p interactions:',scalscp
 422       endif
 423       r0_corr=cutoff_corr-delt_corr
 424       do i=1,20
 425         aad(i,1)=scalscp*aad(i,1)
 426         aad(i,2)=scalscp*aad(i,2)
 427         bad(i,1)=scalscp*bad(i,1)
 428         bad(i,2)=scalscp*bad(i,2)
 429       enddo
 430       call rescale_weights(t_bath)
 431       if(me.eq.king.or..not.out1file)
 432      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 433      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 434      &  wturn4,wturn6
 435    22 format (/'Energy-term weights (scaled):'//
 436      & 'WSCC=   ',f10.6,' (SC-SC)'/
 437      & 'WSCP=   ',f10.6,' (SC-p)'/
 438      & 'WELEC=  ',f10.6,' (p-p electr)'/
 439      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 440      & 'WBOND=  ',f10.6,' (stretching)'/
 441      & 'WANG=   ',f10.6,' (bending)'/
 442      & 'WSCLOC= ',f10.6,' (SC local)'/
 443      & 'WTOR=   ',f10.6,' (torsional)'/
 444      & 'WTORD=  ',f10.6,' (double torsional)'/
 445      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 446      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 447      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 448      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 449      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 450      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 451      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 452      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 453      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 454       if(me.eq.king.or..not.out1file)
 455      & write (iout,*) "Reference temperature for weights calculation:",
 456      &  temp0
 457       call reada(weightcard,"D0CM",d0cm,3.78d0)
 458       call reada(weightcard,"AKCM",akcm,15.1d0)
 459       call reada(weightcard,"AKTH",akth,11.0d0)
 460       call reada(weightcard,"AKCT",akct,12.0d0)
 461       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 462       call reada(weightcard,"V2SS",v2ss,7.61d0)
 463       call reada(weightcard,"V3SS",v3ss,13.7d0)
 464       call reada(weightcard,"EBR",ebr,-5.50D0)
 465       if(me.eq.king.or..not.out1file) then
 466        write (iout,*) "Parameters of the SS-bond potential:"
 467        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 468      & " AKCT",akct
 469        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 470        write (iout,*) "EBR",ebr
 471        print *,'indpdb=',indpdb,' pdbref=',pdbref
 472       endif
 473       if (indpdb.gt.0 .or. pdbref) then
 474         read(inp,'(a)') pdbname(1)
 475         call pdbinput(1)
 476         call contact(.false.,ncont_ref,icont_ref,co)
 477       endif
 478       if (indpdb.eq.0) then
 479 C Read sequence if not taken from the pdb file.
 480         read (inp,*) nres
 481 c        print *,'nres=',nres
 482         if (iscode.gt.0) then
 483           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 484         else
 485           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 486         endif
 487 C Convert sequence to numeric code
 488         do i=1,nres
 489           itype(i)=rescode(i,sequence(i),iscode)
 490         enddo
 491 C Assign initial virtual bond lengths
 492         do i=2,nres
 493           vbld(i)=vbl
 494           vbld_inv(i)=vblinv
 495         enddo
 496         do i=2,nres-1
 497           vbld(i+nres)=dsc(itype(i))
 498           vbld_inv(i+nres)=dsc_inv(itype(i))
 499 c          write (iout,*) "i",i," itype",itype(i),
 500 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 501         enddo
 502       endif
 503 c      print *,nres
 504 c      print '(20i4)',(itype(i),i=1,nres)
 505       do i=1,nres
 506 #ifdef PROCOR
 507         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
 508 #else
 509         if (itype(i).eq.21) then
 510 #endif
 511           itel(i)=0
 512 #ifdef PROCOR
 513         else if (itype(i+1).ne.20) then
 514 #else
 515         else if (itype(i).ne.20) then
 516 #endif
 517           itel(i)=1
 518         else
 519           itel(i)=2
 520         endif
 521       enddo
 522       if(me.eq.king.or..not.out1file)then
 523        write (iout,*) "ITEL"
 524        do i=1,nres-1
 525          write (iout,*) i,itype(i),itel(i)
 526        enddo
 527        print *,'Call Read_Bridge.'
 528       endif
 529       call read_bridge
 530 C 8/13/98 Set limits to generating the dihedral angles
 531       do i=1,nres
 532         phibound(1,i)=-pi
 533         phibound(2,i)=pi
 534       enddo
 535       read (inp,*) ndih_constr
 536       if (ndih_constr.gt.0) then
 537         read (inp,*) ftors
 538         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 539         if(me.eq.king.or..not.out1file)then
 540          write (iout,*)
 541      &   'There are',ndih_constr,' constraints on phi angles.'
 542          do i=1,ndih_constr
 543           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 544          enddo
 545         endif
 546         do i=1,ndih_constr
 547           phi0(i)=deg2rad*phi0(i)
 548           drange(i)=deg2rad*drange(i)
 549         enddo
 550         if(me.eq.king.or..not.out1file)
 551      &   write (iout,*) 'FTORS',ftors
 552         do i=1,ndih_constr
 553           ii = idih_constr(i)
 554           phibound(1,ii) = phi0(i)-drange(i)
 555           phibound(2,ii) = phi0(i)+drange(i)
 556         enddo
 557       endif
 558       nnt=1
 559 #ifdef MPI
 560       if (me.eq.king) then
 561 #endif
 562        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 563        do i=1,nres
 564          write (iout,'(a3,i5,2f10.1)')
 565      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 566        enddo
 567 #ifdef MP
 568       endif
 569 #endif
 570       nct=nres
 571 cd      print *,'NNT=',NNT,' NCT=',NCT
 572       if (itype(1).eq.21) nnt=2
 573       if (itype(nres).eq.21) nct=nct-1
 574       if (pdbref) then
 575         if(me.eq.king.or..not.out1file)
 576      &   write (iout,'(a,i3)') 'nsup=',nsup
 577         nstart_seq=nnt
 578         if (nsup.le.(nct-nnt+1)) then
 579           do i=0,nct-nnt+1-nsup
 580             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 581               nstart_seq=nnt+i
 582               goto 111
 583             endif
 584           enddo
 585           write (iout,'(a)')
 586      &            'Error - sequences to be superposed do not match.'
 587           stop
 588         else
 589           do i=0,nsup-(nct-nnt+1)
 590             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 591      &      then
 592               nstart_sup=nstart_sup+i
 593               nsup=nct-nnt+1
 594               goto 111
 595             endif
 596           enddo
 597           write (iout,'(a)')
 598      &            'Error - sequences to be superposed do not match.'
 599         endif
 600   111   continue
 601         if (nsup.eq.0) nsup=nct-nnt
 602         if (nstart_sup.eq.0) nstart_sup=nnt
 603         if (nstart_seq.eq.0) nstart_seq=nnt
 604         if(me.eq.king.or..not.out1file)
 605      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 606      &                 ' nstart_seq=',nstart_seq
 607       endif
 608 c--- Zscore rms -------
 609       if (nz_start.eq.0) nz_start=nnt
 610       if (nz_end.eq.0 .and. nsup.gt.0) then
 611         nz_end=nnt+nsup-1
 612       else if (nz_end.eq.0) then
 613         nz_end=nct
 614       endif
 615       if(me.eq.king.or..not.out1file)then
 616        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 617        write (iout,*) 'IZ_SC=',iz_sc
 618       endif
 619 c----------------------
 620       call init_int_table
 621       if (refstr) then
 622         if (.not.pdbref) then
 623           call read_angles(inp,*38)
 624           goto 39
 625    38     write (iout,'(a)') 'Error reading reference structure.'
 626 #ifdef MPI
 627           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 628           stop 'Error reading reference structure'
 629 #endif
 630    39     call chainbuild
 631           call setup_var
 632 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 633           nstart_sup=nnt
 634           nstart_seq=nnt
 635           nsup=nct-nnt+1
 636           do i=1,2*nres
 637             do j=1,3
 638               cref(j,i)=c(j,i)
 639             enddo
 640           enddo
 641           call contact(.true.,ncont_ref,icont_ref,co)
 642         endif
 643         if(me.eq.king.or..not.out1file)
 644      &   write (iout,*) 'Contact order:',co
 645         if (pdbref) then
 646         if(me.eq.king.or..not.out1file)
 647      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 648         do i=1,ncont_ref
 649           do j=1,2
 650             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 651           enddo
 652           if(me.eq.king.or..not.out1file)
 653      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 654      &     icont_ref(1,i),' ',
 655      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 656         enddo
 657         endif
 658       endif
 659 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 660       if (constr_dist.gt.0) then
 661         call read_dist_constr
 662       endif
 663       if (nhpb.gt.0) call hpb_partition
 664 c      write (iout,*) "After read_dist_constr nhpb",nhpb
 665 c      call flush(iout)
 666       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 667      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 668      &    modecalc.ne.10) then
 669 C If input structure hasn't been supplied from the PDB file read or generate
 670 C initial geometry.
 671         if (iranconf.eq.0 .and. .not. extconf) then
 672           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 673      &     write (iout,'(a)') 'Initial geometry will be read in.'
 674           if (read_cart) then
 675             read(inp,'(8f10.5)',end=36,err=36)
 676      &       ((c(l,k),l=1,3),k=1,nres),
 677      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 678             call int_from_cart1(.false.)
 679             do i=1,nres-1
 680               do j=1,3
 681                 dc(j,i)=c(j,i+1)-c(j,i)
 682                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 683               enddo
 684             enddo
 685             do i=nnt,nct
 686               if (itype(i).ne.10) then
 687                 do j=1,3
 688                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 689                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 690                 enddo
 691               endif
 692             enddo
 693             return
 694           else
 695             call read_angles(inp,*36)
 696           endif
 697           goto 37
 698    36     write (iout,'(a)') 'Error reading angle file.'
 699 #ifdef MPI
 700           call mpi_finalize( MPI_COMM_WORLD,IERR )
 701 #endif
 702           stop 'Error reading angle file.'
 703    37     continue
 704         else if (extconf) then
 705          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 706      &    write (iout,'(a)') 'Extended chain initial geometry.'
 707          do i=3,nres
 708           theta(i)=90d0*deg2rad
 709          enddo
 710          do i=4,nres
 711           phi(i)=180d0*deg2rad
 712          enddo
 713          do i=2,nres-1
 714           alph(i)=110d0*deg2rad
 715          enddo
 716          do i=2,nres-1
 717           omeg(i)=-120d0*deg2rad
 718          enddo
 719         else
 720           if(me.eq.king.or..not.out1file)
 721      &     write (iout,'(a)') 'Random-generated initial geometry.'
 722
 723
 724 #ifdef MPI
 725           if (me.eq.king  .or. fg_rank.eq.0 .and. (
 726      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
 727 #endif
 728             do itrial=1,100
 729               itmp=1
 730               call gen_rand_conf(itmp,*30)
 731               goto 40
 732    30         write (iout,*) 'Failed to generate random conformation',
 733      &          ', itrial=',itrial
 734               write (*,*) 'Processor:',me,
 735      &          ' Failed to generate random conformation',
 736      &          ' itrial=',itrial
 737               call intout
 738
 739 #ifdef AIX
 740               call flush_(iout)
 741 #else
 742               call flush(iout)
 743 #endif
 744             enddo
 745             write (iout,'(a,i3,a)') 'Processor:',me,
 746      &        ' error in generating random conformation.'
 747             write (*,'(a,i3,a)') 'Processor:',me,
 748      &        ' error in generating random conformation.'
 749             call flush(iout)
 750 #ifdef MPI
 751             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
 752    40       continue
 753           endif
 754 #else
 755    40     continue
 756 #endif
 757         endif
 758       elseif (modecalc.eq.4) then
 759         read (inp,'(a)') intinname
 760         open (intin,file=intinname,status='old',err=333)
 761         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
 762      &  write (iout,'(a)') 'intinname',intinname
 763         write (*,'(a)') 'Processor',myrank,' intinname',intinname
 764         goto 334
 765   333   write (iout,'(2a)') 'Error opening angle file ',intinname
 766 #ifdef MPI
 767         call MPI_Finalize(MPI_COMM_WORLD,IERR)
 768 #endif
 769         stop 'Error opening angle file.'
 770   334   continue
 771
 772       endif
 773 C Generate distance constraints, if the PDB structure is to be regularized.
 774       if (nthread.gt.0) then
 775         call read_threadbase
 776       endif
 777       call setup_var
 778       if (me.eq.king .or. .not. out1file)
 779      & call intout
 780       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
 781         write (iout,'(/a,i3,a)')
 782      &  'The chain contains',ns,' disulfide-bridging cysteines.'
 783         write (iout,'(20i4)') (iss(i),i=1,ns)
 784         write (iout,'(/a/)') 'Pre-formed links are:'
 785         do i=1,nss
 786           i1=ihpb(i)-nres
 787           i2=jhpb(i)-nres
 788           it1=itype(i1)
 789           it2=itype(i2)
 790           if (me.eq.king.or..not.out1file)
 791      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 792      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 793      &    ebr,forcon(i)
 794         enddo
 795         write (iout,'(a)')
 796       endif
 797       if (i2ndstr.gt.0) call secstrp2dihc
 798 c      call geom_to_var(nvar,x)
 799 c      call etotal(energia(0))
 800 c      call enerprint(energia(0))
 801 c      call briefout(0,etot)
 802 c      stop
 803 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
 804 cd    write (iout,'(a)') 'Variable list:'
 805 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
 806       if (indpdb.gt.0) then
 807         do i=1,10000
 808           read(1,'(a)',end=555,err=555) pdbname(i)
 809         enddo
 810   555   npdbfile=i-1
 811         write (iout,*) "PDB files"
 812         do i=1,npdbfile
 813           write (iout,'(i6,a)') i,pdbname(i)(:ilen(pdbname(i)))
 814         enddo
 815       endif
 816 #ifdef MPI
 817       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
 818      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
 819      &  'Processor',myrank,': end reading molecular data.'
 820 #endif
 821       return
 822       end
 823 c--------------------------------------------------------------------------
 824       logical function seq_comp(itypea,itypeb,length)
 825       implicit none
 826       integer length,itypea(length),itypeb(length)
 827       integer i
 828       do i=1,length
 829         if (itypea(i).ne.itypeb(i)) then
 830           seq_comp=.false.
 831           return
 832         endif
 833       enddo
 834       seq_comp=.true.
 835       return
 836       end
 837 c-----------------------------------------------------------------------------
 838       subroutine read_bridge
 839 C Read information about disulfide bridges.
 840       implicit real*8 (a-h,o-z)
 841       include 'DIMENSIONS'
 842 #ifdef MPI
 843       include 'mpif.h'
 844 #endif
 845       include 'COMMON.IOUNITS'
 846       include 'COMMON.GEO'
 847       include 'COMMON.VAR'
 848       include 'COMMON.INTERACT'
 849       include 'COMMON.LOCAL'
 850       include 'COMMON.NAMES'
 851       include 'COMMON.CHAIN'
 852       include 'COMMON.FFIELD'
 853       include 'COMMON.SBRIDGE'
 854       include 'COMMON.HEADER'
 855       include 'COMMON.CONTROL'
 856       include 'COMMON.DBASE'
 857       include 'COMMON.THREAD'
 858       include 'COMMON.TIME1'
 859       include 'COMMON.SETUP'
 860 C Read bridging residues.
 861       read (inp,*) ns,(iss(i),i=1,ns)
 862       print *,'ns=',ns
 863       if(me.eq.king.or..not.out1file)
 864      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
 865 C Check whether the specified bridging residues are cystines.
 866       do i=1,ns
 867         if (itype(iss(i)).ne.1) then
 868           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
 869      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 870      &   ' can form a disulfide bridge?!!!'
 871           write (*,'(2a,i3,a)')
 872      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
 873      &   ' can form a disulfide bridge?!!!'
 874 #ifdef MPI
 875          call MPI_Finalize(MPI_COMM_WORLD,ierror)
 876          stop
 877 #endif
 878         endif
 879       enddo
 880 C Read preformed bridges.
 881       if (ns.gt.0) then
 882       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 883       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
 884       if (nss.gt.0) then
 885         nhpb=nss
 886 C Check if the residues involved in bridges are in the specified list of
 887 C bridging residues.
 888         do i=1,nss
 889           do j=1,i-1
 890             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 891      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 892               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 893      &      ' contains residues present in other pairs.'
 894               write (*,'(a,i3,a)') 'Disulfide pair',i,
 895      &      ' contains residues present in other pairs.'
 896 #ifdef MPI
 897               call MPI_Finalize(MPI_COMM_WORLD,ierror)
 898               stop
 899 #endif
 900             endif
 901           enddo
 902           do j=1,ns
 903             if (ihpb(i).eq.iss(j)) goto 10
 904           enddo
 905           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 906    10     continue
 907           do j=1,ns
 908             if (jhpb(i).eq.iss(j)) goto 20
 909           enddo
 910           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 911    20     continue
 912           dhpb(i)=dbr
 913           forcon(i)=fbr
 914         enddo
 915         do i=1,nss
 916           ihpb(i)=ihpb(i)+nres
 917           jhpb(i)=jhpb(i)+nres
 918         enddo
 919       endif
 920       endif
 921       return
 922       end
 923 c----------------------------------------------------------------------------
 924       subroutine read_x(kanal,*)
 925       implicit real*8 (a-h,o-z)
 926       include 'DIMENSIONS'
 927       include 'COMMON.GEO'
 928       include 'COMMON.VAR'
 929       include 'COMMON.CHAIN'
 930       include 'COMMON.IOUNITS'
 931       include 'COMMON.CONTROL'
 932       include 'COMMON.LOCAL'
 933       include 'COMMON.INTERACT'
 934 c Read coordinates from input
 935 c
 936       read(kanal,'(8f10.5)',end=10,err=10)
 937      &  ((c(l,k),l=1,3),k=1,nres),
 938      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
 939       do j=1,3
 940         c(j,nres+1)=c(j,1)
 941         c(j,2*nres)=c(j,nres)
 942       enddo
 943       call int_from_cart1(.false.)
 944       do i=1,nres-1
 945         do j=1,3
 946           dc(j,i)=c(j,i+1)-c(j,i)
 947           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
 948         enddo
 949       enddo
 950       do i=nnt,nct
 951         if (itype(i).ne.10) then
 952           do j=1,3
 953             dc(j,i+nres)=c(j,i+nres)-c(j,i)
 954             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
 955           enddo
 956         endif
 957       enddo
 958
 959       return
 960    10 return1
 961       end
 962 c----------------------------------------------------------------------------
 963       subroutine read_threadbase
 964       implicit real*8 (a-h,o-z)
 965       include 'DIMENSIONS'
 966       include 'COMMON.IOUNITS'
 967       include 'COMMON.GEO'
 968       include 'COMMON.VAR'
 969       include 'COMMON.INTERACT'
 970       include 'COMMON.LOCAL'
 971       include 'COMMON.NAMES'
 972       include 'COMMON.CHAIN'
 973       include 'COMMON.FFIELD'
 974       include 'COMMON.SBRIDGE'
 975       include 'COMMON.HEADER'
 976       include 'COMMON.CONTROL'
 977       include 'COMMON.DBASE'
 978       include 'COMMON.THREAD'
 979       include 'COMMON.TIME1'
 980 C Read pattern database for threading.
 981       read (icbase,*) nseq
 982       do i=1,nseq
 983         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
 984      &   nres_base(2,i),nres_base(3,i)
 985         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
 986      &   nres_base(1,i))
 987 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
 988 c    &   nres_base(2,i),nres_base(3,i)
 989 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
 990 c    &   nres_base(1,i))
 991       enddo
 992       close (icbase)
 993       if (weidis.eq.0.0D0) weidis=0.1D0
 994       do i=nnt,nct
 995         do j=i+2,nct
 996           nhpb=nhpb+1
 997           ihpb(nhpb)=i
 998           jhpb(nhpb)=j
 999           forcon(nhpb)=weidis
1000         enddo
1001       enddo
1002       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1003       write (iout,'(a,i5)') 'nexcl: ',nexcl
1004       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1005       return
1006       end
1007 c------------------------------------------------------------------------------
1008       subroutine setup_var
1009       implicit real*8 (a-h,o-z)
1010       include 'DIMENSIONS'
1011       include 'COMMON.IOUNITS'
1012       include 'COMMON.GEO'
1013       include 'COMMON.VAR'
1014       include 'COMMON.INTERACT'
1015       include 'COMMON.LOCAL'
1016       include 'COMMON.NAMES'
1017       include 'COMMON.CHAIN'
1018       include 'COMMON.FFIELD'
1019       include 'COMMON.SBRIDGE'
1020       include 'COMMON.HEADER'
1021       include 'COMMON.CONTROL'
1022       include 'COMMON.DBASE'
1023       include 'COMMON.THREAD'
1024       include 'COMMON.TIME1'
1025 C Set up variable list.
1026       ntheta=nres-2
1027       nphi=nres-3
1028       nvar=ntheta+nphi
1029       nside=0
1030       do i=2,nres-1
1031         if (itype(i).ne.10) then
1032           nside=nside+1
1033           ialph(i,1)=nvar+nside
1034           ialph(nside,2)=i
1035         endif
1036       enddo
1037       if (indphi.gt.0) then
1038         nvar=nphi
1039       else if (indback.gt.0) then
1040         nvar=nphi+ntheta
1041       else
1042         nvar=nvar+2*nside
1043       endif
1044 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1045       return
1046       end
1047 c----------------------------------------------------------------------------
1048       subroutine gen_dist_constr
1049 C Generate CA distance constraints.
1050       implicit real*8 (a-h,o-z)
1051       include 'DIMENSIONS'
1052       include 'COMMON.IOUNITS'
1053       include 'COMMON.GEO'
1054       include 'COMMON.VAR'
1055       include 'COMMON.INTERACT'
1056       include 'COMMON.LOCAL'
1057       include 'COMMON.NAMES'
1058       include 'COMMON.CHAIN'
1059       include 'COMMON.FFIELD'
1060       include 'COMMON.SBRIDGE'
1061       include 'COMMON.HEADER'
1062       include 'COMMON.CONTROL'
1063       include 'COMMON.DBASE'
1064       include 'COMMON.THREAD'
1065       include 'COMMON.TIME1'
1066       dimension itype_pdb(maxres)
1067       common /pizda/ itype_pdb
1068       character*2 iden
1069 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1070 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1071 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1072 cd     & ' nsup',nsup
1073       do i=nstart_sup,nstart_sup+nsup-1
1074 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1075 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1076         do j=i+2,nstart_sup+nsup-1
1077           nhpb=nhpb+1
1078           ihpb(nhpb)=i+nstart_seq-nstart_sup
1079           jhpb(nhpb)=j+nstart_seq-nstart_sup
1080           forcon(nhpb)=weidis
1081           dhpb(nhpb)=dist(i,j)
1082         enddo
1083       enddo
1084 cd      write (iout,'(a)') 'Distance constraints:'
1085 cd      do i=nss+1,nhpb
1086 cd        ii=ihpb(i)
1087 cd        jj=jhpb(i)
1088 cd        iden='CA'
1089 cd        if (ii.gt.nres) then
1090 cd          iden='SC'
1091 cd          ii=ii-nres
1092 cd          jj=jj-nres
1093 cd        endif
1094 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1095 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1096 cd     &  dhpb(i),forcon(i)
1097 cd      enddo
1098       return
1099       end
1100 c----------------------------------------------------------------------------
1101       subroutine map_read
1102       implicit real*8 (a-h,o-z)
1103       include 'DIMENSIONS'
1104       include 'COMMON.MAP'
1105       include 'COMMON.IOUNITS'
1106       character*3 angid(4) /'THE','PHI','ALP','OME'/
1107       character*80 mapcard,ucase
1108       do imap=1,nmap
1109         read (inp,'(a)') mapcard
1110         mapcard=ucase(mapcard)
1111         if (index(mapcard,'PHI').gt.0) then
1112           kang(imap)=1
1113         else if (index(mapcard,'THE').gt.0) then
1114           kang(imap)=2
1115         else if (index(mapcard,'ALP').gt.0) then
1116           kang(imap)=3
1117         else if (index(mapcard,'OME').gt.0) then
1118           kang(imap)=4
1119         else
1120           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1121           stop 'Error - illegal variable spec in MAP card.'
1122         endif
1123         call readi (mapcard,'RES1',res1(imap),0)
1124         call readi (mapcard,'RES2',res2(imap),0)
1125         if (res1(imap).eq.0) then
1126           res1(imap)=res2(imap)
1127         else if (res2(imap).eq.0) then
1128           res2(imap)=res1(imap)
1129         endif
1130         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1131           write (iout,'(a)')
1132      &    'Error - illegal definition of variable group in MAP.'
1133           stop 'Error - illegal definition of variable group in MAP.'
1134         endif
1135         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1136         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1137         call readi(mapcard,'NSTEP',nstep(imap),0)
1138         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1139           write (iout,'(a)')
1140      &     'Illegal boundary and/or step size specification in MAP.'
1141           stop 'Illegal boundary and/or step size specification in MAP.'
1142         endif
1143       enddo ! imap
1144       return
1145       end
1146 c----------------------------------------------------------------------------
1147       subroutine read_minim
1148       implicit real*8 (a-h,o-z)
1149       include 'DIMENSIONS'
1150       include 'COMMON.MINIM'
1151       include 'COMMON.IOUNITS'
1152       character*80 ucase
1153       character*320 minimcard
1154       call card_concat(minimcard)
1155       call readi(minimcard,'MAXMIN',maxmin,2000)
1156       call readi(minimcard,'MAXFUN',maxfun,5000)
1157       call readi(minimcard,'MINMIN',minmin,maxmin)
1158       call readi(minimcard,'MINFUN',minfun,maxmin)
1159       call reada(minimcard,'TOLF',tolf,1.0D-2)
1160       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1161       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1162       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1163       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1164       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1165      &         'Options in energy minimization:'
1166       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1167      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1168      & 'MinMin:',MinMin,' MinFun:',MinFun,
1169      & ' TolF:',TolF,' RTolF:',RTolF
1170       return
1171       end
1172 c----------------------------------------------------------------------------
1173       subroutine read_angles(kanal,*)
1174       implicit real*8 (a-h,o-z)
1175       include 'DIMENSIONS'
1176       include 'COMMON.GEO'
1177       include 'COMMON.VAR'
1178       include 'COMMON.CHAIN'
1179       include 'COMMON.IOUNITS'
1180       include 'COMMON.CONTROL'
1181 c Read angles from input
1182 c
1183        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1184        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1185        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1186        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1187
1188        do i=1,nres
1189 c 9/7/01 avoid 180 deg valence angle
1190         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1191 c
1192         theta(i)=deg2rad*theta(i)
1193         phi(i)=deg2rad*phi(i)
1194         alph(i)=deg2rad*alph(i)
1195         omeg(i)=deg2rad*omeg(i)
1196        enddo
1197       return
1198    10 return1
1199       end
1200 c----------------------------------------------------------------------------
1201       subroutine reada(rekord,lancuch,wartosc,default)
1202       implicit none
1203       character*(*) rekord,lancuch
1204       double precision wartosc,default
1205       integer ilen,iread
1206       external ilen
1207       iread=index(rekord,lancuch)
1208       if (iread.eq.0) then
1209         wartosc=default
1210         return
1211       endif
1212       iread=iread+ilen(lancuch)+1
1213       read (rekord(iread:),*,err=10,end=10) wartosc
1214       return
1215   10  wartosc=default
1216       return
1217       end
1218 c----------------------------------------------------------------------------
1219       subroutine readi(rekord,lancuch,wartosc,default)
1220       implicit none
1221       character*(*) rekord,lancuch
1222       integer wartosc,default
1223       integer ilen,iread
1224       external ilen
1225       iread=index(rekord,lancuch)
1226       if (iread.eq.0) then
1227         wartosc=default
1228         return
1229       endif
1230       iread=iread+ilen(lancuch)+1
1231       read (rekord(iread:),*,err=10,end=10) wartosc
1232       return
1233   10  wartosc=default
1234       return
1235       end
1236 c----------------------------------------------------------------------------
1237       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1238       implicit none
1239       integer dim,i
1240       integer tablica(dim),default
1241       character*(*) rekord,lancuch
1242       character*80 aux
1243       integer ilen,iread
1244       external ilen
1245       do i=1,dim
1246         tablica(i)=default
1247       enddo
1248       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1249       if (iread.eq.0) return
1250       iread=iread+ilen(lancuch)+1
1251       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1252    10 return
1253       end
1254 c----------------------------------------------------------------------------
1255       subroutine multreada(rekord,lancuch,tablica,dim,default)
1256       implicit none
1257       integer dim,i
1258       double precision tablica(dim),default
1259       character*(*) rekord,lancuch
1260       character*80 aux
1261       integer ilen,iread
1262       external ilen
1263       do i=1,dim
1264         tablica(i)=default
1265       enddo
1266       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1267       if (iread.eq.0) return
1268       iread=iread+ilen(lancuch)+1
1269       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1270    10 return
1271       end
1272 c----------------------------------------------------------------------------
1273       subroutine card_concat(card)
1274       implicit real*8 (a-h,o-z)
1275       include 'DIMENSIONS'
1276       include 'COMMON.IOUNITS'
1277       character*(*) card
1278       character*80 karta,ucase
1279       external ilen
1280       read (inp,'(a)') karta
1281       karta=ucase(karta)
1282       card=' '
1283       do while (karta(80:80).eq.'&')
1284         card=card(:ilen(card)+1)//karta(:79)
1285         read (inp,'(a)') karta
1286         karta=ucase(karta)
1287       enddo
1288       card=card(:ilen(card)+1)//karta
1289       return
1290       end
1291 c---------------------------------------------------------------------------------
1292       subroutine read_fragments
1293       implicit real*8 (a-h,o-z)
1294       include 'DIMENSIONS'
1295 #ifdef MPI
1296       include 'mpif.h'
1297 #endif
1298       include 'COMMON.SETUP'
1299       include 'COMMON.CHAIN'
1300       include 'COMMON.IOUNITS'
1301       include 'COMMON.MD'
1302       include 'COMMON.CONTROL'
1303       read(inp,*) nset,nfrag,npair,nfrag_back
1304       if(me.eq.king.or..not.out1file)
1305      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
1306      &  " nfrag_back",nfrag_back
1307       do iset=1,nset
1308          read(inp,*) mset(iset)
1309        do i=1,nfrag
1310          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
1311      &     qinfrag(i,iset)
1312          if(me.eq.king.or..not.out1file)
1313      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
1314      &     ifrag(2,i,iset), qinfrag(i,iset)
1315        enddo
1316        do i=1,npair
1317         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
1318      &    qinpair(i,iset)
1319         if(me.eq.king.or..not.out1file)
1320      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
1321      &    ipair(2,i,iset), qinpair(i,iset)
1322        enddo
1323        do i=1,nfrag_back
1324         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
1325      &     wfrag_back(3,i,iset),
1326      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
1327         if(me.eq.king.or..not.out1file)
1328      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
1329      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
1330        enddo
1331       enddo
1332       return
1333       end
1334 c-------------------------------------------------------------------------------
1335       subroutine read_dist_constr
1336       implicit real*8 (a-h,o-z)
1337       include 'DIMENSIONS'
1338 #ifdef MPI
1339       include 'mpif.h'
1340 #endif
1341       include 'COMMON.SETUP'
1342       include 'COMMON.CONTROL'
1343       include 'COMMON.CHAIN'
1344       include 'COMMON.IOUNITS'
1345       include 'COMMON.SBRIDGE'
1346       integer ifrag_(2,100),ipair_(2,100)
1347       double precision wfrag_(100),wpair_(100)
1348       character*500 controlcard
1349 c      write (iout,*) "Calling read_dist_constr"
1350 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
1351 c      call flush(iout)
1352       call card_concat(controlcard)
1353       call readi(controlcard,"NFRAG",nfrag_,0)
1354       call readi(controlcard,"NPAIR",npair_,0)
1355       call readi(controlcard,"NDIST",ndist_,0)
1356       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
1357       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
1358       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
1359       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
1360       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
1361 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
1362 c      write (iout,*) "IFRAG"
1363 c      do i=1,nfrag_
1364 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1365 c      enddo
1366 c      write (iout,*) "IPAIR"
1367 c      do i=1,npair_
1368 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
1369 c      enddo
1370       if (.not.refstr .and. nfrag.gt.0) then
1371         write (iout,*)
1372      &  "ERROR: no reference structure to compute distance restraints"
1373         write (iout,*)
1374      &  "Restraints must be specified explicitly (NDIST=number)"
1375         stop
1376       endif
1377       if (nfrag.lt.2 .and. npair.gt.0) then
1378         write (iout,*) "ERROR: Less than 2 fragments specified",
1379      &   " but distance restraints between pairs requested"
1380         stop
1381       endif
1382       call flush(iout)
1383       do i=1,nfrag_
1384         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
1385         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
1386      &    ifrag_(2,i)=nstart_sup+nsup-1
1387 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
1388         call flush(iout)
1389         if (wfrag_(i).gt.0.0d0) then
1390         do j=ifrag_(1,i),ifrag_(2,i)-1
1391           do k=j+1,ifrag_(2,i)
1392             write (iout,*) "j",j," k",k
1393             ddjk=dist(j,k)
1394             if (constr_dist.eq.1) then
1395             nhpb=nhpb+1
1396             ihpb(nhpb)=j
1397             jhpb(nhpb)=k
1398               dhpb(nhpb)=ddjk
1399             forcon(nhpb)=wfrag_(i)
1400             else if (constr_dist.eq.2) then
1401               if (ddjk.le.dist_cut) then
1402                 nhpb=nhpb+1
1403                 ihpb(nhpb)=j
1404                 jhpb(nhpb)=k
1405                 dhpb(nhpb)=ddjk
1406                 forcon(nhpb)=wfrag_(i)
1407               endif
1408             else
1409               nhpb=nhpb+1
1410               ihpb(nhpb)=j
1411               jhpb(nhpb)=k
1412               dhpb(nhpb)=ddjk
1413               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
1414             endif
1415 #ifdef MPI
1416             if (.not.out1file .or. me.eq.king)
1417      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1418      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1419 #else
1420             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
1421      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1422 #endif
1423           enddo
1424         enddo
1425         endif
1426       enddo
1427       do i=1,npair_
1428         if (wpair_(i).gt.0.0d0) then
1429         ii = ipair_(1,i)
1430         jj = ipair_(2,i)
1431         if (ii.gt.jj) then
1432           itemp=ii
1433           ii=jj
1434           jj=itemp
1435         endif
1436         do j=ifrag_(1,ii),ifrag_(2,ii)
1437           do k=ifrag_(1,jj),ifrag_(2,jj)
1438             nhpb=nhpb+1
1439             ihpb(nhpb)=j
1440             jhpb(nhpb)=k
1441             forcon(nhpb)=wpair_(i)
1442             dhpb(nhpb)=dist(j,k)
1443 #ifdef MPI
1444             if (.not.out1file .or. me.eq.king)
1445      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1446      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1447 #else
1448             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
1449      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
1450 #endif
1451           enddo
1452         enddo
1453         endif
1454       enddo
1455       do i=1,ndist_
1456         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
1457      &     ibecarb(i),forcon(nhpb+1)
1458         if (forcon(nhpb+1).gt.0.0d0) then
1459           nhpb=nhpb+1
1460           if (ibecarb(i).gt.0) then
1461             ihpb(i)=ihpb(i)+nres
1462             jhpb(i)=jhpb(i)+nres
1463           endif
1464           if (dhpb(nhpb).eq.0.0d0)
1465      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
1466         endif
1467       enddo
1468 #ifdef MPI
1469       if (.not.out1file .or. me.eq.king) then
1470 #endif
1471       do i=1,nhpb
1472           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
1473      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
1474       enddo
1475       call flush(iout)
1476 #ifdef MPI
1477       endif
1478 #endif
1479       return
1480       end
1481 c-------------------------------------------------------------------------------
1482 #ifdef WINIFL
1483       subroutine flush(iu)
1484       return
1485       end
1486 #endif
1487 #ifdef AIX
1488       subroutine flush(iu)
1489       call flush_(iu)
1490       return
1491       end
1492 #endif
1493 c------------------------------------------------------------------------------
1494       subroutine copy_to_tmp(source)
1495       include "DIMENSIONS"
1496       include "COMMON.IOUNITS"
1497       character*(*) source
1498       character* 256 tmpfile
1499       integer ilen
1500       external ilen
1501       logical ex
1502       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
1503       inquire(file=tmpfile,exist=ex)
1504       if (ex) then
1505         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
1506      &   " to temporary directory..."
1507         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
1508         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
1509       endif
1510       return
1511       end
1512 c------------------------------------------------------------------------------
1513       subroutine move_from_tmp(source)
1514       include "DIMENSIONS"
1515       include "COMMON.IOUNITS"
1516       character*(*) source
1517       integer ilen
1518       external ilen
1519       write (*,*) "Moving ",source(:ilen(source)),
1520      & " from temporary directory to working directory"
1521       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
1522       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
1523       return
1524       end
1525 c------------------------------------------------------------------------------
1526       subroutine random_init(seed)
1527 C
1528 C Initialize random number generator
1529 C
1530       implicit real*8 (a-h,o-z)
1531       include 'DIMENSIONS'
1532 #ifdef AMD64
1533       integer*8 iseedi8
1534 #endif
1535 #ifdef MPI
1536       include 'mpif.h'
1537       logical OKRandom, prng_restart
1538       real*8  r1
1539       integer iseed_array(4)
1540 #endif
1541       include 'COMMON.IOUNITS'
1542       include 'COMMON.TIME1'
1543       include 'COMMON.THREAD'
1544       include 'COMMON.SBRIDGE'
1545       include 'COMMON.CONTROL'
1546       include 'COMMON.MCM'
1547       include 'COMMON.MAP'
1548       include 'COMMON.HEADER'
1549 csa      include 'COMMON.CSA'
1550       include 'COMMON.CHAIN'
1551       include 'COMMON.MUCA'
1552       include 'COMMON.MD'
1553       include 'COMMON.FFIELD'
1554       include 'COMMON.SETUP'
1555       iseed=-dint(dabs(seed))
1556       if (iseed.eq.0) then
1557         write (iout,'(/80(1h*)/20x,a/80(1h*))')
1558      &    'Random seed undefined. The program will stop.'
1559         write (*,'(/80(1h*)/20x,a/80(1h*))')
1560      &    'Random seed undefined. The program will stop.'
1561 #ifdef MPI
1562         call mpi_finalize(mpi_comm_world,ierr)
1563 #endif
1564         stop 'Bad random seed.'
1565       endif
1566 #ifdef MPI
1567       if (fg_rank.eq.0) then
1568       seed=seed*(me+1)+1
1569 #ifdef AMD64
1570       iseedi8=dint(seed)
1571       if(me.eq.king .or. .not. out1file)
1572      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1573       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
1574       OKRandom = prng_restart(me,iseedi8)
1575 #else
1576       do i=1,4
1577        tmp=65536.0d0**(4-i)
1578        iseed_array(i) = dint(seed/tmp)
1579        seed=seed-iseed_array(i)*tmp
1580       enddo
1581       if(me.eq.king .or. .not. out1file)
1582      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
1583      &                 (iseed_array(i),i=1,4)
1584       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
1585      &                 (iseed_array(i),i=1,4)
1586       OKRandom = prng_restart(me,iseed_array)
1587 #endif
1588       if (OKRandom) then
1589         r1=ran_number(0.0D0,1.0D0)
1590         if(me.eq.king .or. .not. out1file)
1591      &   write (iout,*) 'ran_num',r1
1592         if (r1.lt.0.0d0) OKRandom=.false.
1593       endif
1594       if (.not.OKRandom) then
1595         write (iout,*) 'PRNG IS NOT WORKING!!!'
1596         print *,'PRNG IS NOT WORKING!!!'
1597         if (me.eq.0) then
1598          call flush(iout)
1599          call mpi_abort(mpi_comm_world,error_msg,ierr)
1600          stop
1601         else
1602          write (iout,*) 'too many processors for parallel prng'
1603          write (*,*) 'too many processors for parallel prng'
1604          call flush(iout)
1605          stop
1606         endif
1607       endif
1608       endif
1609 #else
1610       call vrndst(iseed)
1611       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
1612 #endif
1613       return
1614       end
1615 c--------------------------------------------------------------------------
1616       subroutine pdbinput(iconf)
1617       implicit real*8 (a-h,o-z)
1618       include 'DIMENSIONS'
1619 #ifdef MPI
1620       include 'mpif.h'
1621       integer error_msg
1622 #endif
1623       include 'COMMON.IOUNITS'
1624       include 'COMMON.GEO'
1625       include 'COMMON.VAR'
1626       include 'COMMON.INTERACT'
1627       include 'COMMON.LOCAL'
1628       include 'COMMON.NAMES'
1629       include 'COMMON.CHAIN'
1630       include 'COMMON.FFIELD'
1631       include 'COMMON.SBRIDGE'
1632       include 'COMMON.HEADER'
1633       include 'COMMON.CONTROL'
1634       include 'COMMON.DBASE'
1635       include 'COMMON.THREAD'
1636       include 'COMMON.CONTACTS'
1637       include 'COMMON.TORCNSTR'
1638       include 'COMMON.TIME1'
1639       include 'COMMON.BOUNDS'
1640       include 'COMMON.MD'
1641       include 'COMMON.REMD'
1642       include 'COMMON.SETUP'
1643       character*4 sequence(maxres)
1644       integer rescode
1645       double precision x(maxvar)
1646       dimension itype_pdb(maxres)
1647       common /pizda/ itype_pdb
1648       logical seq_comp,fail
1649       double precision energia(0:n_ene)
1650       integer ilen
1651       external ilen
1652       if(me.eq.king.or..not.out1file)
1653      & write (iout,'(2a)') 'PDB data will be read from file ',
1654      & pdbname(iconf)(:ilen(pdbname(iconf)))
1655       open(ipdbin,file=pdbname(iconf),status='old',err=33)
1656       goto 34
1657   33  write (iout,'(a)') 'Error opening PDB file.'
1658       stop
1659   34  continue
1660 c      print *,'Begin reading pdb data'
1661       call readpdb
1662 c      print *,'Finished reading pdb data'
1663       if(me.eq.king.or..not.out1file)
1664      & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
1665      & ' nstart_sup=',nstart_sup
1666       do i=1,nres
1667         itype_pdb(i)=itype(i)
1668       enddo
1669       close (ipdbin)
1670       nnt=nstart_sup
1671       nct=nstart_sup+nsup-1
1672
1673       if (sideadd) then
1674 C Following 2 lines for diagnostics; comment out if not needed
1675         write (iout,*) "Before sideadd"
1676         call intout
1677         if (me.eq.king.or..not.out1file)
1678      &    write(iout,*)'Adding sidechains'
1679         maxsi=1000
1680         do i=2,nres-1
1681           iti=itype(i)
1682           if (iti.ne.10) then
1683             nsi=0
1684             fail=.true.
1685             do while (fail.and.nsi.le.maxsi)
1686               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
1687               nsi=nsi+1
1688             enddo
1689             if(fail) write(iout,*)'Adding sidechain failed for res ',
1690      &            i,' after ',nsi,' trials'
1691           endif
1692         enddo
1693 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
1694         call chainbuild
1695 C Following 2 lines for diagnostics; comment out if not needed
1696         write (iout,*) "After sideadd"
1697         call intout
1698       endif
1699       return
1700       end