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