source/cluster/wham/src-M/readrtns.F

   1       subroutine read_control
   2 C
   3 C Read molecular data
   4 C
   5       implicit none
   6       include 'DIMENSIONS'
   7       include 'sizesclu.dat'
   8       include 'COMMON.IOUNITS'
   9       include 'COMMON.TIME1'
  10       include 'COMMON.SBRIDGE'
  11       include 'COMMON.CONTROL'
  12       include 'COMMON.CLUSTER'
  13       include 'COMMON.CHAIN'
  14       include 'COMMON.HEADER'
  15       include 'COMMON.FFIELD'
  16       include 'COMMON.FREE'
  17       include 'COMMON.INTERACT'
  18       include "COMMON.SPLITELE"
  19       character*320 controlcard,ucase
  20 #ifdef MPL
  21       include 'COMMON.INFO'
  22 #endif
  23       integer i,i1,i2,it1,it2
  24
  25       read (INP,'(a80)') titel
  26       call card_concat(controlcard)
  27
  28       call readi(controlcard,'NRES',nres,0)
  29       call readi(controlcard,'RESCALE',rescale_mode,2)
  30       call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0)
  31       write (iout,*) "DISTCHAINMAX",distchainmax
  32 C Reading the dimensions of box in x,y,z coordinates
  33       call reada(controlcard,'BOXX',boxxsize,100.0d0)
  34       call reada(controlcard,'BOXY',boxysize,100.0d0)
  35       call reada(controlcard,'BOXZ',boxzsize,100.0d0)
  36 c Cutoff range for interactions
  37       call reada(controlcard,"R_CUT",r_cut,15.0d0)
  38       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
  39       call readi(controlcard,'PDBOUT',outpdb,0)
  40       call readi(controlcard,'MOL2OUT',outmol2,0)
  41       refstr=(index(controlcard,'REFSTR').gt.0)
  42       write (iout,*) "REFSTR",refstr
  43       pdbref=(index(controlcard,'PDBREF').gt.0)
  44       iscode=index(controlcard,'ONE_LETTER')
  45       tree=(index(controlcard,'MAKE_TREE').gt.0)
  46       with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
  47       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  48       write (iout,*) "with_dihed_constr ",with_dihed_constr,
  49      & " CONSTR_DIST",constr_dist
  50       with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
  51       write (iout,*) "with_theta_constr ",with_theta_constr
  52       call flush(iout)
  53       min_var=(index(controlcard,'MINVAR').gt.0)
  54       plot_tree=(index(controlcard,'PLOT_TREE').gt.0)
  55       punch_dist=(index(controlcard,'PUNCH_DIST').gt.0)
  56       call readi(controlcard,'NCUT',ncut,1)
  57       call readi(controlcard,'SYM',symetr,1)
  58       write (iout,*) 'sym', symetr
  59       call readi(controlcard,'NSTART',nstart,0)
  60       call readi(controlcard,'NEND',nend,0)
  61       call reada(controlcard,'ECUT',ecut,10.0d0)
  62       call reada(controlcard,'PROB',prob_limit,0.99d0)
  63       write (iout,*) "Probability limit",prob_limit
  64       lgrp=(index(controlcard,'LGRP').gt.0)
  65       caonly=(index(controlcard,'CA_ONLY').gt.0)
  66       print_dist=(index(controlcard,'PRINT_DIST').gt.0)
  67       call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0)
  68       call readi(controlcard,'IOPT',iopt,2)
  69       lside = index(controlcard,"SIDE").gt.0
  70       efree = index(controlcard,"EFREE").gt.0
  71       call readi(controlcard,'NTEMP',nT,1)
  72       write (iout,*) "nT",nT
  73       call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0)
  74       write (iout,*) "nT",nT
  75       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  76       do i=1,nT
  77         beta_h(i)=1.0d0/(1.987D-3*beta_h(i))
  78       enddo
  79       write (iout,*) 'beta_h',(beta_h(i),i=1,nT)
  80       lprint_cart=index(controlcard,"PRINT_CART") .gt.0
  81       lprint_int=index(controlcard,"PRINT_INT") .gt.0
  82       if (min_var) iopt=1
  83       return
  84       end
  85 c--------------------------------------------------------------------------
  86       subroutine molread
  87 C
  88 C Read molecular data.
  89 C
  90       implicit none
  91       include 'DIMENSIONS'
  92       include 'COMMON.IOUNITS'
  93       include 'COMMON.GEO'
  94       include 'COMMON.VAR'
  95       include 'COMMON.INTERACT'
  96       include 'COMMON.LOCAL'
  97       include 'COMMON.NAMES'
  98       include 'COMMON.CHAIN'
  99       include 'COMMON.FFIELD'
 100       include 'COMMON.SBRIDGE'
 101       include 'COMMON.HEADER'
 102       include 'COMMON.CONTROL'
 103       include 'COMMON.CONTACTS'
 104       include 'COMMON.TIME1'
 105       include 'COMMON.TORCNSTR'
 106 #ifdef MPL
 107       include 'COMMON.INFO'
 108 #endif
 109       character*4 sequence(maxres)
 110       character*800 weightcard
 111       integer rescode
 112       double precision x(maxvar)
 113       integer itype_pdb(maxres)
 114       logical seq_comp
 115       integer i,j,kkk,i1,i2,it1,it2
 116 C
 117 C Body
 118 C
 119 C Read weights of the subsequent energy terms.
 120       call card_concat(weightcard)
 121       write(iout,*) weightcard
 122       call reada(weightcard,'WSC',wsc,1.0d0)
 123       call reada(weightcard,'WLONG',wsc,wsc)
 124       write(iout,*) "WLONG=",wsc
 125       call reada(weightcard,'WSCP',wscp,1.0d0)
 126       call reada(weightcard,'WELEC',welec,1.0D0)
 127       call reada(weightcard,'WVDWPP',wvdwpp,welec)
 128       call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 129       call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 130       call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 131       call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 132       call reada(weightcard,'WTURN3',wturn3,1.0D0)
 133       call reada(weightcard,'WTURN4',wturn4,1.0D0)
 134       call reada(weightcard,'WTURN6',wturn6,1.0D0)
 135       call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 136       call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 137       call reada(weightcard,'WBOND',wbond,1.0D0)
 138       call reada(weightcard,'WTOR',wtor,1.0D0)
 139       call reada(weightcard,'WTORD',wtor_d,1.0D0)
 140       call reada(weightcard,'WANG',wang,1.0D0)
 141       call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 142       call reada(weightcard,'SCAL14',scal14,0.4D0)
 143       call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 144       call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 145       call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 146       if (index(weightcard,'SOFT').gt.0) ipot=6
 147       call reada(weightcard,"D0CM",d0cm,3.78d0)
 148       call reada(weightcard,"AKCM",akcm,15.1d0)
 149       call reada(weightcard,"AKTH",akth,11.0d0)
 150       call reada(weightcard,"AKCT",akct,12.0d0)
 151       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 152       call reada(weightcard,"V2SS",v2ss,7.61d0)
 153       call reada(weightcard,"V3SS",v3ss,13.7d0)
 154       call reada(weightcard,"EBR",ebr,-5.50D0)
 155       write (iout,*) "atriss",atriss
 156       call reada(weightcard,"ATRISS",atriss,0.301D0)
 157       call reada(weightcard,"BTRISS",btriss,0.021D0)
 158       call reada(weightcard,"CTRISS",ctriss,1.001D0)
 159       call reada(weightcard,"DTRISS",dtriss,1.001D0)
 160       write(iout,*) "after"
 161       write (iout,*) "ATRISS=", atriss
 162       write (iout,*) "BTRISS=", btriss
 163       write (iout,*) "CTRISS=", ctriss
 164       write (iout,*) "DTRISS=", dtriss
 165       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 166       do i=1,maxres
 167         dyn_ss_mask(i)=.false.
 168       enddo
 169       do i=1,maxres-1
 170         do j=i+1,maxres
 171           dyn_ssbond_ij(i,j)=1.0d300
 172         enddo
 173       enddo
 174       call reada(weightcard,"HT",Ht,0.0D0)
 175       if (dyn_ss) then
 176         ss_depth=ebr/wsc-0.25*eps(1,1)
 177         Ht=Ht/wsc-0.25*eps(1,1)
 178         akcm=akcm*wstrain/wsc
 179         akth=akth*wstrain/wsc
 180         akct=akct*wstrain/wsc
 181         v1ss=v1ss*wstrain/wsc
 182         v2ss=v2ss*wstrain/wsc
 183         v3ss=v3ss*wstrain/wsc
 184       else
 185         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 186       endif
 187       write (iout,'(/a)') "Disulfide bridge parameters:"
 188       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
 189       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
 190       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
 191       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
 192      & ' v3ss:',v3ss
 193
 194 C 12/1/95 Added weight for the multi-body term WCORR
 195       call reada(weightcard,'WCORRH',wcorr,1.0D0)
 196       if (wcorr4.gt.0.0d0) wcorr=wcorr4
 197       weights(1)=wsc
 198       weights(2)=wscp
 199       weights(3)=welec
 200       weights(4)=wcorr
 201       weights(5)=wcorr5
 202       weights(6)=wcorr6
 203       weights(7)=wel_loc
 204       weights(8)=wturn3
 205       weights(9)=wturn4
 206       weights(10)=wturn6
 207       weights(11)=wang
 208       weights(12)=wscloc
 209       weights(13)=wtor
 210       weights(14)=wtor_d
 211       weights(15)=wstrain
 212       weights(16)=wvdwpp
 213       weights(17)=wbond
 214       weights(18)=scal14
 215       weights(19)=wsccor
 216       write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 217      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wturn3,
 218      &  wturn4,wturn6,wsccor
 219    10 format (/'Energy-term weights (unscaled):'//
 220      & 'WSC=    ',f10.6,' (SC-SC)'/
 221      & 'WSCP=   ',f10.6,' (SC-p)'/
 222      & 'WELEC=  ',f10.6,' (p-p electr)'/
 223      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 224      & 'WBOND=  ',f10.6,' (stretching)'/
 225      & 'WANG=   ',f10.6,' (bending)'/
 226      & 'WSCLOC= ',f10.6,' (SC local)'/
 227      & 'WTOR=   ',f10.6,' (torsional)'/
 228      & 'WTORD=  ',f10.6,' (double torsional)'/
 229      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 230      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 231      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 232      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 233      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 234      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 235      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 236      & 'WTURN6= ',f10.6,' (turns, 6th order)'/
 237      & 'WSCCOR= ',f10.6,' (SC-backbone torsinal correalations)')
 238
 239       if (wcorr4.gt.0.0d0) then
 240         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 241      &   'between contact pairs of peptide groups'
 242         write (iout,'(2(a,f5.3/))')
 243      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 244      &  'Range of quenching the correlation terms:',2*delt_corr
 245       else if (wcorr.gt.0.0d0) then
 246         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 247      &   'between contact pairs of peptide groups'
 248       endif
 249       write (iout,'(a,f8.3)')
 250      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 251       write (iout,'(a,f8.3)')
 252      &  'General scaling factor of SC-p interactions:',scalscp
 253       r0_corr=cutoff_corr-delt_corr
 254       do i=1,20
 255         aad(i,1)=scalscp*aad(i,1)
 256         aad(i,2)=scalscp*aad(i,2)
 257         bad(i,1)=scalscp*bad(i,1)
 258         bad(i,2)=scalscp*bad(i,2)
 259       enddo
 260
 261       call flush(iout)
 262       print *,'indpdb=',indpdb,' pdbref=',pdbref
 263
 264 C Read sequence if not taken from the pdb file.
 265       if (iscode.gt.0) then
 266         read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 267       else
 268         read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 269       endif
 270 C Convert sequence to numeric code
 271       do i=1,nres
 272         itype(i)=rescode(i,sequence(i),iscode)
 273       enddo
 274       print *,nres
 275       print '(20i4)',(itype(i),i=1,nres)
 276
 277       do i=1,nres
 278 #ifdef PROCOR
 279         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 280 #else
 281         if (itype(i).eq.ntyp1) then
 282 #endif
 283           itel(i)=0
 284 #ifdef PROCOR
 285         else if (iabs(itype(i+1)).ne.20) then
 286 #else
 287         else if (iabs(itype(i)).ne.20) then
 288 #endif
 289           itel(i)=1
 290         else
 291           itel(i)=2
 292         endif
 293       enddo
 294       write (iout,*) "ITEL"
 295       do i=1,nres-1
 296         write (iout,*) i,itype(i),itel(i)
 297       enddo
 298
 299       print *,'Call Read_Bridge.'
 300       call read_bridge
 301 C this fragment reads diheadral constrains
 302       if (with_dihed_constr) then
 303
 304       read (inp,*) ndih_constr
 305       if (ndih_constr.gt.0) then
 306 C        read (inp,*) ftors
 307 C        write (iout,*) 'FTORS',ftors
 308 C ftors is the force constant for torsional quartic constrains
 309         read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i),
 310      &                i=1,ndih_constr)
 311         write (iout,*)
 312      &   'There are',ndih_constr,' constraints on phi angles.'
 313         do i=1,ndih_constr
 314           write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i),
 315      &  ftors(i)
 316         enddo
 317         do i=1,ndih_constr
 318           phi0(i)=deg2rad*phi0(i)
 319           drange(i)=deg2rad*drange(i)
 320         enddo
 321       endif ! endif ndif_constr.gt.0
 322       endif ! with_dihed_constr
 323       if (with_theta_constr) then
 324 C with_theta_constr is keyword allowing for occurance of theta constrains
 325       read (inp,*) ntheta_constr
 326 C ntheta_constr is the number of theta constrains
 327       if (ntheta_constr.gt.0) then
 328 C        read (inp,*) ftors
 329         read (inp,*) (itheta_constr(i),theta_constr0(i),
 330      &  theta_drange(i),for_thet_constr(i),
 331      &  i=1,ntheta_constr)
 332 C the above code reads from 1 to ntheta_constr
 333 C itheta_constr(i) residue i for which is theta_constr
 334 C theta_constr0 the global minimum value
 335 C theta_drange is range for which there is no energy penalty
 336 C for_thet_constr is the force constant for quartic energy penalty
 337 C E=k*x**4
 338 C        if(me.eq.king.or..not.out1file)then
 339          write (iout,*)
 340      &   'There are',ntheta_constr,' constraints on phi angles.'
 341          do i=1,ntheta_constr
 342           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 343      &    theta_drange(i),
 344      &    for_thet_constr(i)
 345          enddo
 346 C        endif
 347         do i=1,ntheta_constr
 348           theta_constr0(i)=deg2rad*theta_constr0(i)
 349           theta_drange(i)=deg2rad*theta_drange(i)
 350         enddo
 351 C        if(me.eq.king.or..not.out1file)
 352 C     &   write (iout,*) 'FTORS',ftors
 353 C        do i=1,ntheta_constr
 354 C          ii = itheta_constr(i)
 355 C          thetabound(1,ii) = phi0(i)-drange(i)
 356 C          thetabound(2,ii) = phi0(i)+drange(i)
 357 C        enddo
 358       endif ! ntheta_constr.gt.0
 359       endif! with_theta_constr
 360
 361       nnt=1
 362       nct=nres
 363       print *,'NNT=',NNT,' NCT=',NCT
 364       if (itype(1).eq.ntyp1) nnt=2
 365       if (itype(nres).eq.ntyp1) nct=nct-1
 366       if (nstart.lt.nnt) nstart=nnt
 367       if (nend.gt.nct .or. nend.eq.0) nend=nct
 368       write (iout,*) "nstart",nstart," nend",nend
 369       nres0=nres
 370 c      if (pdbref) then
 371 c        read(inp,'(a)') pdbfile
 372 c        write (iout,'(2a)') 'PDB data will be read from file ',pdbfile
 373 c        open(ipdbin,file=pdbfile,status='old',err=33)
 374 c        goto 34
 375 c  33    write (iout,'(a)') 'Error opening PDB file.'
 376 c        stop
 377 c  34    continue
 378 c        print *,'Begin reading pdb data'
 379 c        call readpdb
 380 c        print *,'Finished reading pdb data'
 381 c        write (iout,'(a,i3,a,i3)')'nsup=',nsup,' nstart_sup=',nstart_sup
 382 c        do i=1,nres
 383 c          itype_pdb(i)=itype(i)
 384 c        enddo
 385 c        close (ipdbin)
 386 c        write (iout,'(a,i3)') 'nsup=',nsup
 387 c        nstart_seq=nnt
 388 c        if (nsup.le.(nct-nnt+1)) then
 389 c          do i=0,nct-nnt+1-nsup
 390 c            if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 391 c              nstart_seq=nnt+i
 392 c              goto 111
 393 c            endif
 394 c          enddo
 395 c          write (iout,'(a)')
 396 c     &            'Error - sequences to be superposed do not match.'
 397 c          stop
 398 c        else
 399 c          do i=0,nsup-(nct-nnt+1)
 400 c            if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 401 c     &      then
 402 c              nstart_sup=nstart_sup+i
 403 c              nsup=nct-nnt+1
 404 c              goto 111
 405 c            endif
 406 c          enddo
 407 c          write (iout,'(a)')
 408 c     &            'Error - sequences to be superposed do not match.'
 409 c        endif
 410 c  111   continue
 411 c        write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 412 c     &                 ' nstart_seq=',nstart_seq
 413 c      endif
 414       call init_int_table
 415       call setup_var
 416       write (iout,*) "molread: REFSTR",refstr
 417       if (refstr) then
 418         if (.not.pdbref) then
 419           call read_angles(inp,*38)
 420           goto 39
 421    38     write (iout,'(a)') 'Error reading reference structure.'
 422 #ifdef MPL
 423           call mp_stopall(Error_Msg)
 424 #else
 425           stop 'Error reading reference structure'
 426 #endif
 427    39     call chainbuild
 428           nstart_sup=nnt
 429           nstart_seq=nnt
 430           nsup=nct-nnt+1
 431           kkk=1
 432           do i=1,2*nres
 433             do j=1,3
 434               cref(j,i,kkk)=c(j,i)
 435             enddo
 436           enddo
 437         endif
 438         call contact(.true.,ncont_ref,icont_ref)
 439       endif
 440        if (ns.gt.0) then
 441 C        write (iout,'(/a,i3,a)')
 442 C     &  'The chain contains',ns,' disulfide-bridging cysteines.'
 443         write (iout,'(20i4)') (iss(i),i=1,ns)
 444        if (dyn_ss) then
 445           write(iout,*)"Running with dynamic disulfide-bond formation"
 446        else
 447         write (iout,'(/a/)') 'Pre-formed links are:'
 448         do i=1,nss
 449           i1=ihpb(i)-nres
 450           i2=jhpb(i)-nres
 451           it1=itype(i1)
 452           it2=itype(i2)
 453           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
 454      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
 455      &    ebr,forcon(i)
 456         enddo
 457         write (iout,'(a)')
 458        endif
 459       endif
 460       if (ns.gt.0.and.dyn_ss) then
 461           do i=nss+1,nhpb
 462             ihpb(i-nss)=ihpb(i)
 463             jhpb(i-nss)=jhpb(i)
 464             forcon(i-nss)=forcon(i)
 465             dhpb(i-nss)=dhpb(i)
 466           enddo
 467           nhpb=nhpb-nss
 468           nss=0
 469           call hpb_partition
 470           do i=1,ns
 471             dyn_ss_mask(iss(i))=.true.
 472           enddo
 473       endif
 474 c Read distance restraints
 475       if (constr_dist.gt.0) then
 476         call read_dist_constr
 477         call hpb_partition
 478       endif
 479       return
 480       end
 481 c-----------------------------------------------------------------------------
 482       logical function seq_comp(itypea,itypeb,length)
 483       implicit none
 484       integer length,itypea(length),itypeb(length)
 485       integer i
 486       do i=1,length
 487         if (itypea(i).ne.itypeb(i)) then
 488           seq_comp=.false.
 489           return
 490         endif
 491       enddo
 492       seq_comp=.true.
 493       return
 494       end
 495 c-----------------------------------------------------------------------------
 496       subroutine read_bridge
 497 C Read information about disulfide bridges.
 498       implicit none
 499       include 'DIMENSIONS'
 500       include 'COMMON.IOUNITS'
 501       include 'COMMON.GEO'
 502       include 'COMMON.VAR'
 503       include 'COMMON.INTERACT'
 504       include 'COMMON.LOCAL'
 505       include 'COMMON.NAMES'
 506       include 'COMMON.CHAIN'
 507       include 'COMMON.FFIELD'
 508       include 'COMMON.SBRIDGE'
 509       include 'COMMON.HEADER'
 510       include 'COMMON.CONTROL'
 511       include 'COMMON.TIME1'
 512 #ifdef MPL
 513       include 'COMMON.INFO'
 514 #endif
 515       integer i,j
 516 C Read bridging residues.
 517       read (inp,*) ns,(iss(i),i=1,ns)
 518       print *,'ns=',ns
 519 C Check whether the specified bridging residues are cystines.
 520       do i=1,ns
 521         if (itype(iss(i)).ne.1) then
 522           write (iout,'(2a,i3,a)')
 523      &   'Do you REALLY think that the residue ',
 524      &    restyp(itype(iss(i))),i,
 525      &   ' can form a disulfide bridge?!!!'
 526           write (*,'(2a,i3,a)')
 527      &   'Do you REALLY think that the residue ',
 528      &   restyp(itype(iss(i))),i,
 529      &   ' can form a disulfide bridge?!!!'
 530 #ifdef MPL
 531          call mp_stopall(error_msg)
 532 #else
 533          stop
 534 #endif
 535         endif
 536       enddo
 537 C Read preformed bridges.
 538       if (ns.gt.0) then
 539       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
 540       if (nss.gt.0) then
 541         nhpb=nss
 542 C Check if the residues involved in bridges are in the specified list of
 543 C bridging residues.
 544         do i=1,nss
 545           do j=1,i-1
 546             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
 547      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
 548               write (iout,'(a,i3,a)') 'Disulfide pair',i,
 549      &      ' contains residues present in other pairs.'
 550               write (*,'(a,i3,a)') 'Disulfide pair',i,
 551      &      ' contains residues present in other pairs.'
 552 #ifdef MPL
 553               call mp_stopall(error_msg)
 554 #else
 555               stop
 556 #endif
 557             endif
 558           enddo
 559           do j=1,ns
 560             if (ihpb(i).eq.iss(j)) goto 10
 561           enddo
 562           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 563    10     continue
 564           do j=1,ns
 565             if (jhpb(i).eq.iss(j)) goto 20
 566           enddo
 567           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
 568    20     continue
 569 C          dhpb(i)=dbr
 570 C          forcon(i)=fbr
 571         enddo
 572         do i=1,nss
 573           ihpb(i)=ihpb(i)+nres
 574           jhpb(i)=jhpb(i)+nres
 575         enddo
 576       endif
 577       endif
 578       return
 579       end
 580 c----------------------------------------------------------------------------
 581       subroutine read_angles(kanal,*)
 582       implicit none
 583       include 'DIMENSIONS'
 584       include 'COMMON.GEO'
 585       include 'COMMON.VAR'
 586       include 'COMMON.CHAIN'
 587       include 'COMMON.IOUNITS'
 588       integer i,kanal
 589       read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
 590       read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
 591       read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
 592       read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
 593       do i=1,nres
 594         theta(i)=deg2rad*theta(i)
 595         phi(i)=deg2rad*phi(i)
 596         alph(i)=deg2rad*alph(i)
 597         omeg(i)=deg2rad*omeg(i)
 598       enddo
 599       return
 600    10 return1
 601       end
 602 c----------------------------------------------------------------------------
 603       subroutine reada(rekord,lancuch,wartosc,default)
 604       implicit none
 605       character*(*) rekord,lancuch
 606       double precision wartosc,default
 607       integer ilen,iread
 608       external ilen
 609       iread=index(rekord,lancuch)
 610       if (iread.eq.0) then
 611         wartosc=default
 612         return
 613       endif
 614       iread=iread+ilen(lancuch)+1
 615       read (rekord(iread:),*) wartosc
 616       return
 617       end
 618 c----------------------------------------------------------------------------
 619       subroutine multreada(rekord,lancuch,tablica,dim,default)
 620       implicit none
 621       integer dim,i
 622       double precision tablica(dim),default
 623       character*(*) rekord,lancuch
 624       integer ilen,iread
 625       external ilen
 626       do i=1,dim
 627         tablica(i)=default
 628       enddo
 629       iread=index(rekord,lancuch)
 630       if (iread.eq.0) return
 631       iread=iread+ilen(lancuch)+1
 632       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 633    10 return
 634       end
 635 c----------------------------------------------------------------------------
 636       subroutine readi(rekord,lancuch,wartosc,default)
 637       implicit none
 638       character*(*) rekord,lancuch
 639       integer wartosc,default
 640       integer ilen,iread
 641       external ilen
 642       iread=index(rekord,lancuch)
 643       if (iread.eq.0) then
 644         wartosc=default
 645         return
 646       endif
 647       iread=iread+ilen(lancuch)+1
 648       read (rekord(iread:),*) wartosc
 649       return
 650       end
 651 C----------------------------------------------------------------------
 652       subroutine multreadi(rekord,lancuch,tablica,dim,default)
 653       implicit none
 654       integer dim,i
 655       integer tablica(dim),default
 656       character*(*) rekord,lancuch
 657       character*80 aux
 658       integer ilen,iread
 659       external ilen
 660       do i=1,dim
 661         tablica(i)=default
 662       enddo
 663       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
 664       if (iread.eq.0) return
 665       iread=iread+ilen(lancuch)+1
 666       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
 667    10 return
 668       end
 669
 670 c----------------------------------------------------------------------------
 671       subroutine card_concat(card)
 672       include 'DIMENSIONS'
 673       include 'COMMON.IOUNITS'
 674       character*(*) card
 675       character*80 karta,ucase
 676       external ilen
 677       read (inp,'(a)') karta
 678       karta=ucase(karta)
 679       card=' '
 680       do while (karta(80:80).eq.'&')
 681         card=card(:ilen(card)+1)//karta(:79)
 682         read (inp,'(a)') karta
 683         karta=ucase(karta)
 684       enddo
 685       card=card(:ilen(card)+1)//karta
 686       return
 687       end
 688 c----------------------------------------------------------------------------
 689       subroutine openunits
 690       implicit none
 691       include 'DIMENSIONS'
 692 #ifdef MPI
 693       include "mpif.h"
 694       character*3 liczba
 695       include "COMMON.MPI"
 696 #endif
 697       include 'COMMON.IOUNITS'
 698       include 'COMMON.CONTROL'
 699       integer lenpre,lenpot,ilen
 700       external ilen
 701       character*16 cformat,cprint
 702       character*16 ucase
 703       integer lenint,lenout
 704       call getenv('INPUT',prefix)
 705       call getenv('OUTPUT',prefout)
 706       call getenv('INTIN',prefintin)
 707       call getenv('COORD',cformat)
 708       call getenv('PRINTCOOR',cprint)
 709       call getenv('SCRATCHDIR',scratchdir)
 710       from_bx=.true.
 711       from_cx=.false.
 712       if (index(ucase(cformat),'CX').gt.0) then
 713         from_cx=.true.
 714         from_bx=.false.
 715       endif
 716       from_cart=.true.
 717       lenpre=ilen(prefix)
 718       lenout=ilen(prefout)
 719       lenint=ilen(prefintin)
 720 C Get the names and open the input files
 721       open (inp,file=prefix(:ilen(prefix))//'.inp',status='old')
 722 #ifdef MPI
 723       write (liczba,'(bz,i3.3)') me
 724       outname=prefout(:lenout)//'_clust.out_'//liczba
 725 #else
 726       outname=prefout(:lenout)//'_clust.out'
 727 #endif
 728       if (from_bx) then
 729         intinname=prefintin(:lenint)//'.bx'
 730       else if (from_cx) then
 731         intinname=prefintin(:lenint)//'.cx'
 732       else
 733         intinname=prefintin(:lenint)//'.int'
 734       endif
 735       rmsname=prefintin(:lenint)//'.rms'
 736       open (jplot,file=prefout(:ilen(prefout))//'.tex',
 737      &       status='unknown')
 738       open (jrms,file=rmsname,status='unknown')
 739       open(iout,file=outname,status='unknown')
 740 C Get parameter filenames and open the parameter files.
 741       call getenv('BONDPAR',bondname)
 742       open (ibond,file=bondname,status='old')
 743       call getenv('THETPAR',thetname)
 744       open (ithep,file=thetname,status='old')
 745       call getenv('ROTPAR',rotname)
 746       open (irotam,file=rotname,status='old')
 747       call getenv('TORPAR',torname)
 748       open (itorp,file=torname,status='old')
 749       call getenv('TORDPAR',tordname)
 750       open (itordp,file=tordname,status='old')
 751       call getenv('FOURIER',fouriername)
 752       open (ifourier,file=fouriername,status='old')
 753       call getenv('ELEPAR',elename)
 754       open (ielep,file=elename,status='old')
 755       call getenv('SIDEPAR',sidename)
 756       open (isidep,file=sidename,status='old')
 757       call getenv('SIDEP',sidepname)
 758       open (isidep1,file=sidepname,status="old")
 759       call getenv('SCCORPAR',sccorname)
 760       open (isccor,file=sccorname,status="old")
 761 #ifndef OLDSCP
 762 C
 763 C 8/9/01 In the newest version SCp interaction constants are read from a file
 764 C Use -DOLDSCP to use hard-coded constants instead.
 765 C
 766       call getenv('SCPPAR',scpname)
 767       open (iscpp,file=scpname,status='old')
 768 #endif
 769       return
 770       end
 771       subroutine read_dist_constr
 772       implicit real*8 (a-h,o-z)
 773       include 'DIMENSIONS'
 774 #ifdef MPI
 775       include 'mpif.h'
 776 #endif
 777       include 'COMMON.CONTROL'
 778       include 'COMMON.CHAIN'
 779       include 'COMMON.IOUNITS'
 780       include 'COMMON.SBRIDGE'
 781       integer ifrag_(2,100),ipair_(2,100)
 782       double precision wfrag_(100),wpair_(100)
 783       character*500 controlcard
 784       logical lprn /.true./
 785       write (iout,*) "Calling read_dist_constr"
 786 C      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
 787 C      call flush(iout)
 788       write(iout,*) "TU sie wywalam?"
 789       call card_concat(controlcard)
 790       write (iout,*) controlcard
 791       call flush(iout)
 792       call readi(controlcard,"NFRAG",nfrag_,0)
 793       call readi(controlcard,"NPAIR",npair_,0)
 794       call readi(controlcard,"NDIST",ndist_,0)
 795       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
 796       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
 797       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
 798       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
 799       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
 800       write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 801       write (iout,*) "IFRAG"
 802       do i=1,nfrag_
 803         write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 804       enddo
 805       write (iout,*) "IPAIR"
 806       do i=1,npair_
 807         write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
 808       enddo
 809       call flush(iout)
 810       do i=1,nfrag_
 811         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
 812         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
 813      &    ifrag_(2,i)=nstart_sup+nsup-1
 814 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
 815         call flush(iout)
 816         if (wfrag_(i).gt.0.0d0) then
 817         do j=ifrag_(1,i),ifrag_(2,i)-1
 818           do k=j+1,ifrag_(2,i)
 819             write (iout,*) "j",j," k",k
 820             ddjk=dist(j,k)
 821             if (constr_dist.eq.1) then
 822               nhpb=nhpb+1
 823               ihpb(nhpb)=j
 824               jhpb(nhpb)=k
 825               dhpb(nhpb)=ddjk
 826               forcon(nhpb)=wfrag_(i)
 827             else if (constr_dist.eq.2) then
 828               if (ddjk.le.dist_cut) then
 829                 nhpb=nhpb+1
 830                 ihpb(nhpb)=j
 831                 jhpb(nhpb)=k
 832                 dhpb(nhpb)=ddjk
 833                 forcon(nhpb)=wfrag_(i)
 834               endif
 835             else
 836               nhpb=nhpb+1
 837               ihpb(nhpb)=j
 838               jhpb(nhpb)=k
 839               dhpb(nhpb)=ddjk
 840               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
 841             endif
 842             if (lprn)
 843      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
 844      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 845           enddo
 846         enddo
 847         endif
 848       enddo
 849       do i=1,npair_
 850         if (wpair_(i).gt.0.0d0) then
 851         ii = ipair_(1,i)
 852         jj = ipair_(2,i)
 853         if (ii.gt.jj) then
 854           itemp=ii
 855           ii=jj
 856           jj=itemp
 857         endif
 858         do j=ifrag_(1,ii),ifrag_(2,ii)
 859           do k=ifrag_(1,jj),ifrag_(2,jj)
 860             nhpb=nhpb+1
 861             ihpb(nhpb)=j
 862             jhpb(nhpb)=k
 863             forcon(nhpb)=wpair_(i)
 864             dhpb(nhpb)=dist(j,k)
 865             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 866      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 867           enddo
 868         enddo
 869         endif
 870       enddo
 871       do i=1,ndist_
 872         if (constr_dist.eq.11) then
 873         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
 874      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
 875         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
 876 C        write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 877 C     &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 878         else
 879         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
 880         endif
 881         if (forcon(nhpb+1).gt.0.0d0) then
 882           nhpb=nhpb+1
 883           if (ibecarb(i).gt.0) then
 884             ihpb(i)=ihpb(i)+nres
 885             jhpb(i)=jhpb(i)+nres
 886           endif
 887           if (dhpb(nhpb).eq.0.0d0)
 888      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 889 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 890           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
 891      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
 892         endif
 893 C      endif
 894       enddo
 895       call hpb_partition
 896       call flush(iout)
 897       return
 898       end