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