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