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