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