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