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