2 !-----------------------------------------------------------------------
8 !-----------------------------------------------------------------------------
11 !-----------------------------------------------------------------------------
13 !-----------------------------------------------------------------------------
15 !-----------------------------------------------------------------------------
16 subroutine write_csa_pdb(var,ene,nft,ik,iw_pdb)
19 use geometry_data, only:nres,nvar
20 use geometry, only:var_to_geom,chainbuild
21 use compare, only:secondary2
22 ! implicit real*8 (a-h,o-z)
23 ! include 'DIMENSIONS'
24 ! include 'COMMON.CSA'
25 ! include 'COMMON.BANK'
26 ! include 'COMMON.VAR'
27 ! include 'COMMON.IOUNITS'
28 ! include 'COMMON.MINIM'
29 ! include 'COMMON.SETUP'
30 ! include 'COMMON.GEO'
31 ! include 'COMMON.CHAIN'
32 ! include 'COMMON.LOCAL'
33 ! include 'COMMON.INTERACT'
34 ! include 'COMMON.NAMES'
35 ! include 'COMMON.SBRIDGE'
36 integer :: lenpre,lenpot !,ilen
38 real(kind=8),dimension(nvar) :: var !(maxvar) (maxvar=6*maxres)
39 character(len=50) :: titelloc
40 character(len=3) :: zahl
41 real(kind=8),dimension(mxch*(mxch+1)/2+1) :: ene
43 integer :: nft,ik,iw_pdb
46 if(ene(1).lt.eglob_csa) then
49 call numstr(nglob_csa,zahl)
51 call var_to_geom(nvar,var)
53 call secondary2(.false.)
56 open(icsa_pdb,file=prefix(:lenpre)//'@'//zahl//'.pdb')
59 write(titelloc,'(a2,i3,a3,i9,a3,i6)') &
60 'GM',nglob_csa,' e ',nft,' m ',nmin_csa
62 write(titelloc,'(a2,i3,a3,i9,a3,i6,a5,f5.2,a5,f5.1)') &
63 'GM',nglob_csa,' e ',nft,' m ',nmin_csa,' rms ',&
64 rmsn(ik),' %NC ',pncn(ik)*100
66 call pdbout(eglob_csa,titelloc,icsa_pdb)
71 end subroutine write_csa_pdb
72 !-----------------------------------------------------------------------------
74 !-----------------------------------------------------------------------------
75 subroutine from_pdb(n,idum)
76 ! This subroutine stores the UNRES int variables generated from
77 ! subroutine readpdb into the 1st conformation of in dihang_in.
78 ! Subsequent n-1 conformations of dihang_in have identical values
79 ! of theta and phi as the 1st conformation but random values for
81 ! The array cref (also generated from subroutine readpdb) is stored
82 ! to crefjlee to be used for rmsd calculation in CSA, if necessary.
86 use random, only: ran1
87 ! implicit real*8 (a-h,o-z)
88 ! include 'DIMENSIONS'
89 ! include 'COMMON.IOUNITS'
90 ! include 'COMMON.CHAIN'
91 ! include 'COMMON.VAR'
92 ! include 'COMMON.BANK'
93 ! include 'COMMON.GEO'
95 integer :: n,idum,m,i,j,k,kk,kkk
100 dihang_in(1,j,1,m)=theta(j+1)
101 dihang_in(2,j,1,m)=phi(j+2)
102 dihang_in(3,j,1,m)=alph(j)
103 dihang_in(4,j,1,m)=omeg(j)
105 dihang_in(2,nres-1,1,k)=0.0d0
109 dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
110 dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
111 if(dabs(dihang_in(3,k,1,1)).gt.1.d-6) then
112 dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
113 dihang_in(3,k,1,m)=dihang_in(3,k,1,m)*deg2rad
115 if(dabs(dihang_in(4,k,1,1)).gt.1.d-6) then
116 dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
117 dihang_in(4,k,1,m)=dihang_in(4,k,1,m)*deg2rad
122 ! Store cref to crefjlee (they are in COMMON.CHAIN).
126 crefjlee(kk,k)=cref(kk,k,kkk)
130 open(icsa_native_int,file=csa_native_int,status="old")
132 write(icsa_native_int,*) m,e
133 write(icsa_native_int,200) &
134 (dihang_in(1,k,1,m)*rad2deg,k=2,nres-1)
135 write(icsa_native_int,200) &
136 (dihang_in(2,k,1,m)*rad2deg,k=2,nres-2)
137 write(icsa_native_int,200) &
138 (dihang_in(3,k,1,m)*rad2deg,k=2,nres-1)
139 write(icsa_native_int,200) &
140 (dihang_in(4,k,1,m)*rad2deg,k=2,nres-1)
144 write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
146 close(icsa_native_int)
151 end subroutine from_pdb
152 !-----------------------------------------------------------------------------
153 subroutine from_int(n,mm,idum)
158 use geometry, only:chainbuild,gen_side
159 use energy, only:etotal
161 ! implicit real*8 (a-h,o-z)
162 ! include 'DIMENSIONS'
163 ! include 'COMMON.IOUNITS'
164 ! include 'COMMON.CHAIN'
165 ! include 'COMMON.VAR'
166 ! include 'COMMON.INTERACT'
167 ! include 'COMMON.BANK'
168 ! include 'COMMON.GEO'
169 ! include 'COMMON.CONTACTS'
173 real(kind=8),dimension(0:n_ene) :: energia
175 integer :: n,mm,idum,i,ii,j,m,k,kk,maxcount_fail,icount_fail,maxsi
178 open(icsa_native_int,file=csa_native_int,status="old")
179 read (icsa_native_int,*)
180 call read_angles(icsa_native_int,*10)
182 10 write (iout,'(2a)') "CHUJ NASTAPIL - error in ",&
183 csa_native_int(:ilen(csa_native_int))
187 dihang_in(1,j,1,1)=theta(j+1)
188 dihang_in(2,j,1,1)=phi(j+2)
189 dihang_in(3,j,1,1)=alph(j)
190 dihang_in(4,j,1,1)=omeg(j)
192 dihang_in(2,nres-1,1,1)=0.0d0
194 ! read(icsa_native_int,*) ind,e
195 ! read(icsa_native_int,200) (dihang_in(1,k,1,1),k=2,nres-1)
196 ! read(icsa_native_int,200) (dihang_in(2,k,1,1),k=2,nres-2)
197 ! read(icsa_native_int,200) (dihang_in(3,k,1,1),k=2,nres-1)
198 ! read(icsa_native_int,200) (dihang_in(4,k,1,1),k=2,nres-1)
199 ! dihang_in(2,nres-1,1,1)=0.d0
206 ! dihang_in(1,k,1,m)=dihang_in(1,k,1,1)
207 ! dihang_in(2,k,1,m)=dihang_in(2,k,1,1)
208 ! if(abs(dihang_in(3,k,1,1)).gt.1.d-3) then
209 ! dihang_in(3,k,1,m)=90.d0*ran1(idum)+90.d0
211 ! if(abs(dihang_in(4,k,1,1)).gt.1.d-3) then
212 ! dihang_in(4,k,1,m)=360.d0*ran1(idum)-180.d0
220 DO WHILE (FAIL .AND. ICOUNT_FAIL .LE. MAXCOUNT_FAIL)
223 if (itype(i,1).ne.10) then
224 !d print *,'i=',i,' itype=',itype(i,1),' theta=',theta(i+1)
227 do while (fail .and. ii .le. maxsi)
228 call gen_side(itype(i,1),theta(i+1),alph(i),omeg(i),fail,molnum(i))
235 fail = (energia(0).ge.1.0d20)
236 icount_fail=icount_fail+1
240 if (icount_fail.gt.maxcount_fail) then
242 'Failed to generate non-overlaping near-native conf.',&
247 dihang_in(1,j,1,m)=theta(j+1)
248 dihang_in(2,j,1,m)=phi(j+2)
249 dihang_in(3,j,1,m)=alph(j)
250 dihang_in(4,j,1,m)=omeg(j)
252 dihang_in(2,nres-1,1,m)=0.0d0
256 ! write(icsa_native_int,*) m,e
257 ! write(icsa_native_int,200) (dihang_in(1,k,1,m),k=2,nres-1)
258 ! write(icsa_native_int,200) (dihang_in(2,k,1,m),k=2,nres-2)
259 ! write(icsa_native_int,200) (dihang_in(3,k,1,m),k=2,nres-1)
260 ! write(icsa_native_int,200) (dihang_in(4,k,1,m),k=2,nres-1)
262 ! close(icsa_native_int)
267 ! dihang_in(i,j,1,m)=dihang_in(i,j,1,m)*deg2rad
272 call dihang_to_c(dihang_in(1,1,1,1))
274 ! Store c to cref (they are in COMMON.CHAIN).
277 crefjlee(kk,k)=c(kk,k)
281 call contact(.true.,ncont_ref,icont_ref,co)
284 ! write(icsa_native_int,200) (crefjlee(i,k),i=1,3)
286 close(icsa_native_int)
291 end subroutine from_int
292 !-----------------------------------------------------------------------------
293 subroutine dihang_to_c(aarray)
297 use geometry, only:chainbuild
298 ! implicit real*8 (a-h,o-z)
299 ! include 'DIMENSIONS'
300 ! include 'COMMON.CSA'
301 ! include 'COMMON.BANK'
302 ! include 'COMMON.CHAIN'
303 ! include 'COMMON.GEO'
304 ! include 'COMMON.VAR'
306 real(kind=8),dimension(mxang,nres,mxch) :: aarray !(mxang,maxres,mxch)
309 ! phi(i)=dihang_in(1,i-2,1,1)
312 theta(i+1)=aarray(1,i,1)
313 phi(i+2)=aarray(2,i,1)
314 alph(i)=aarray(3,i,1)
315 omeg(i)=aarray(4,i,1)
321 end subroutine dihang_to_c
322 !-----------------------------------------------------------------------------
324 !-----------------------------------------------------------------------------
326 subroutine cartout(time)
328 subroutine cartoutx(time)
330 use geometry_data, only: c,nres
332 use MD_data, only: potE,t_bath
333 ! implicit real*8 (a-h,o-z)
334 ! include 'DIMENSIONS'
335 ! include 'COMMON.CHAIN'
336 ! include 'COMMON.INTERACT'
337 ! include 'COMMON.NAMES'
338 ! include 'COMMON.IOUNITS'
339 ! include 'COMMON.HEADER'
340 ! include 'COMMON.SBRIDGE'
341 ! include 'COMMON.DISTFIT'
342 ! include 'COMMON.MD'
347 #if defined(AIX) || defined(PGI)
348 open(icart,file=cartname,position="append")
350 open(icart,file=cartname,access="append")
352 write (icart,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
354 write (icart,'(i4,$)') &
355 nss,(idssb(j)+nres,jdssb(j)+nres,j=1,nss)
357 write (icart,'(i4,$)') &
358 nss,(ihpb(j),jhpb(j),j=1,nss)
360 write (icart,'(i4,20f7.4)') nfrag+npair+3*nfrag_back,&
361 (qfrag(i),i=1,nfrag),(qpair(i),i=1,npair),&
362 (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
363 write (icart,'(8f10.5)') &
364 ((c(k,j),k=1,3),j=1,nres),&
365 ((c(k,j+nres),k=1,3),j=nnt,nct)
370 end subroutine cartout
372 end subroutine cartoutx
374 !-----------------------------------------------------------------------------
376 subroutine cartout(time)
377 ! implicit real*8 (a-h,o-z)
378 ! include 'DIMENSIONS'
379 use geometry_data, only: c,nres
381 use MD_data, only: potE,t_bath
385 ! include 'COMMON.SETUP'
387 integer,parameter :: me=0
389 ! include 'COMMON.CHAIN'
390 ! include 'COMMON.INTERACT'
391 ! include 'COMMON.NAMES'
392 ! include 'COMMON.IOUNITS'
393 ! include 'COMMON.HEADER'
394 ! include 'COMMON.SBRIDGE'
395 ! include 'COMMON.DISTFIT'
396 ! include 'COMMON.MD'
400 real(kind=4),dimension(3,2*nres+2) :: xcoord !(3,maxres2+2) (maxres2=2*maxres
405 call xdrfopen_(ixdrf,cartname, "a", iret)
406 call xdrffloat_(ixdrf, real(time), iret)
407 call xdrffloat_(ixdrf, real(potE), iret)
408 call xdrffloat_(ixdrf, real(uconst), iret)
409 call xdrffloat_(ixdrf, real(uconst_back), iret)
410 call xdrffloat_(ixdrf, real(t_bath), iret)
411 call xdrfint_(ixdrf, nss, iret)
414 call xdrfint_(ixdrf, idssb(j)+nres, iret)
415 call xdrfint_(ixdrf, jdssb(j)+nres, iret)
417 call xdrfint_(ixdrf, ihpb(j), iret)
418 call xdrfint_(ixdrf, jhpb(j), iret)
421 call xdrfint_(ixdrf, nfrag+npair+3*nfrag_back, iret)
423 call xdrffloat_(ixdrf, real(qfrag(i)), iret)
426 call xdrffloat_(ixdrf, real(qpair(i)), iret)
429 call xdrffloat_(ixdrf, real(utheta(i)), iret)
430 call xdrffloat_(ixdrf, real(ugamma(i)), iret)
431 call xdrffloat_(ixdrf, real(uscdiff(i)), iret)
434 call xdrfopen(ixdrf,cartname, "a", iret)
435 call xdrffloat(ixdrf, real(time), iret)
436 call xdrffloat(ixdrf, real(potE), iret)
437 call xdrffloat(ixdrf, real(uconst), iret)
438 call xdrffloat(ixdrf, real(uconst_back), iret)
439 call xdrffloat(ixdrf, real(t_bath), iret)
440 call xdrfint(ixdrf, nss, iret)
443 call xdrfint(ixdrf, idssb(j)+nres, iret)
444 call xdrfint(ixdrf, jdssb(j)+nres, iret)
446 call xdrfint(ixdrf, ihpb(j), iret)
447 call xdrfint(ixdrf, jhpb(j), iret)
450 call xdrfint(ixdrf, nfrag+npair+3*nfrag_back, iret)
452 call xdrffloat(ixdrf, real(qfrag(i)), iret)
455 call xdrffloat(ixdrf, real(qpair(i)), iret)
458 call xdrffloat(ixdrf, real(utheta(i)), iret)
459 call xdrffloat(ixdrf, real(ugamma(i)), iret)
460 call xdrffloat(ixdrf, real(uscdiff(i)), iret)
471 xcoord(j,nres+i-nnt+1)=c(j,i+nres)
477 call xdrf3dfcoord_(ixdrf, xcoord, itmp, prec, iret)
478 call xdrfclose_(ixdrf, iret)
480 call xdrf3dfcoord(ixdrf, xcoord, itmp, prec, iret)
481 call xdrfclose(ixdrf, iret)
484 end subroutine cartout
486 !-----------------------------------------------------------------------------
487 subroutine statout(itime)
493 use compare, only:rms_nac_nnc
494 ! implicit real*8 (a-h,o-z)
495 ! include 'DIMENSIONS'
496 ! include 'COMMON.CONTROL'
497 ! include 'COMMON.CHAIN'
498 ! include 'COMMON.INTERACT'
499 ! include 'COMMON.NAMES'
500 ! include 'COMMON.IOUNITS'
501 ! include 'COMMON.HEADER'
502 ! include 'COMMON.SBRIDGE'
503 ! include 'COMMON.DISTFIT'
504 ! include 'COMMON.MD'
505 ! include 'COMMON.REMD'
506 ! include 'COMMON.SETUP'
508 real(kind=8),dimension(0:n_ene) :: energia
509 ! double precision gyrate
511 !el common /gucio/ cm
512 character(len=256) :: line1,line2
513 character(len=4) :: format1,format2
514 character(len=30) :: format
516 integer :: i,ii1,ii2,j
517 real(kind=8) :: rms,frac,frac_nn,co,distance
521 open(istat,file=statname,position="append")
525 open(istat,file=statname,position="append")
527 open(istat,file=statname,access="append")
530 if (AFMlog.gt.0) then
532 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
533 write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
534 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
535 rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
539 !C print *,'A CHUJ',potEcomp(23)
540 write (line1,'(i10,f15.2,7f12.3,i5,$)') &
541 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
542 kinetic_T,t_bath,gyrate(),&
546 else if (selfguide.gt.0) then
549 distance=distance+(c(j,afmend)-c(j,afmbeg))**2
551 distance=dsqrt(distance)
553 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
554 write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
556 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
557 rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
558 distance,potEcomp(23),me
562 !C print *,'A CHUJ',potEcomp(23)
563 write (line1,'(i10,f15.2,8f12.3,i5,$)')&
564 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
565 kinetic_T,t_bath,gyrate(),&
566 distance,potEcomp(23),me
569 else if (velnanoconst.ne.0 ) then
571 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
572 write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
574 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
575 rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
580 !C print *,'A CHUJ',potEcomp(23)
581 write (line1,'(i10,f15.2,8f12.3,i5,$)')&
582 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
583 kinetic_T,t_bath,gyrate(),&
589 call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
590 write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
591 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
592 rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
595 write (line1,'(i10,f15.2,7f12.3,i5,$)') &
596 itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
597 amax,kinetic_T,t_bath,gyrate(),me
601 if(usampl.and.totT.gt.eq_time) then
602 write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back,&
603 (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair),&
604 (utheta(i),ugamma(i),uscdiff(i),i=1,nfrag_back)
605 write(format2,'(a1,i3.3)') "a",23+7*nfrag+7*npair &
611 if (print_compon) then
613 write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
615 write (istat,format) "#","",&
616 (ename(print_order(i)),i=1,nprint_ene)
618 write(format,'(a1,a4,a1,a4,a10)') "(",format1,",",format2,&
620 write (istat,format) line1,line2,&
621 (potEcomp(print_order(i)),i=1,nprint_ene)
623 write(format,'(a1,a4,a1,a4,a1)') "(",format1,",",format2,")"
624 write (istat,format) line1,line2
632 end subroutine statout
633 !-----------------------------------------------------------------------------
635 !-----------------------------------------------------------------------------
641 use muca_md, only:read_muca
642 ! implicit real*8 (a-h,o-z)
643 ! include 'DIMENSIONS'
647 ! include 'COMMON.SETUP'
648 ! include 'COMMON.CONTROL'
649 ! include 'COMMON.SBRIDGE'
650 ! include 'COMMON.IOUNITS'
651 logical :: file_exist
653 ! Read force-field parameters except weights
655 ! Read job setup parameters
657 ! Read force-field parameters except weights
660 ! Read control parameters for energy minimzation if required
661 if (minim) call read_minim
662 ! Read MCM control parameters if required
663 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
664 ! Read MD control parameters if reqjuired
665 if (modecalc.eq.12) call read_MDpar
666 ! Read MREMD control parameters if required
667 if (modecalc.eq.14) then
671 ! Read MUCA control parameters if required
672 if (lmuca) call read_muca
673 ! Read CSA control parameters if required (from fort.40 if exists
674 ! otherwise from general input file)
675 if (modecalc.eq.8) then
676 inquire (file="fort.40",exist=file_exist)
677 if (.not.file_exist) call csaread
679 !fmc if (modecalc.eq.10) call mcmfread
680 ! Read molecule information, molecule geometry, energy-term weights, and
681 ! restraints if requested
683 ! Print restraint information
685 if (.not. out1file .or. me.eq.king) then
688 write (iout,'(a,i5,a)') "The following",nhpb-nss,&
689 " distance constraints have been imposed"
691 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
696 ! print *,"Processor",myrank," leaves READRTNS"
697 ! write(iout,*) "end readrtns"
699 end subroutine readrtns
700 !-----------------------------------------------------------------------------
703 ! Read molecular data.
705 ! use control, only: ilen
711 use MD_data, only: t_bath
713 use compare, only:seq_comp,contact
716 ! implicit real*8 (a-h,o-z)
717 ! include 'DIMENSIONS'
720 integer :: error_msg,ierror,ierr,ierrcode
722 ! include 'COMMON.IOUNITS'
723 ! include 'COMMON.GEO'
724 ! include 'COMMON.VAR'
725 ! include 'COMMON.INTERACT'
726 ! include 'COMMON.LOCAL'
727 ! include 'COMMON.NAMES'
728 ! include 'COMMON.CHAIN'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.SBRIDGE'
731 ! include 'COMMON.HEADER'
732 ! include 'COMMON.CONTROL'
733 ! include 'COMMON.DBASE'
734 ! include 'COMMON.THREAD'
735 ! include 'COMMON.CONTACTS'
736 ! include 'COMMON.TORCNSTR'
737 ! include 'COMMON.TIME1'
738 ! include 'COMMON.BOUNDS'
739 ! include 'COMMON.MD'
740 ! include 'COMMON.SETUP'
741 character(len=4),dimension(:,:),allocatable :: sequence !(maxres,maxmolecules)
743 ! double precision x(maxvar)
744 character(len=256) :: pdbfile
745 character(len=800) :: weightcard
746 character(len=80) :: weightcard_t!,ucase
747 ! integer,dimension(:),allocatable :: itype_pdb !(maxres)
748 ! common /pizda/ itype_pdb
749 logical :: fail !seq_comp,
750 real(kind=8) :: energia(0:n_ene)
754 integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2,mnum
756 real(kind=8),dimension(3,maxres2+2) :: c_alloc
757 real(kind=8),dimension(3,0:maxres2) :: dc_alloc
758 real(kind=8),dimension(:,:), allocatable :: secprob
759 integer,dimension(maxres) :: itype_alloc
761 integer :: iti,nsi,maxsi,itrial,itmp
762 real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,&
763 sigmabet,sumv,nres_temp
764 allocate(weights(n_ene))
765 !-----------------------------
766 allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres
767 allocate(dc(3,0:2*maxres)) !(3,0:maxres2)
768 allocate(itype(maxres,5)) !(maxres)
769 allocate(istype(maxres))
776 !-----------------------------
780 ! Read weights of the subsequent energy terms.
781 call card_concat(weightcard,.true.)
782 call reada(weightcard,'WLONG',wlong,1.0D0)
783 call reada(weightcard,'WSC',wsc,wlong)
784 call reada(weightcard,'WSCP',wscp,wlong)
785 call reada(weightcard,'WELEC',welec,1.0D0)
786 call reada(weightcard,'WVDWPP',wvdwpp,welec)
787 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
788 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
789 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
790 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
791 call reada(weightcard,'WTURN3',wturn3,1.0D0)
792 call reada(weightcard,'WTURN4',wturn4,1.0D0)
793 call reada(weightcard,'WTURN6',wturn6,1.0D0)
794 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
795 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
796 call reada(weightcard,'WVDWPP_NUCL',wvdwpp_nucl,0.0D0)
797 call reada(weightcard,'WELPP',welpp,0.0d0)
798 call reada(weightcard,'WVDWPSB',wvdwpsb,0.0d0)
799 call reada(weightcard,'WELPSB',welpsb,0.0D0)
800 call reada(weightcard,'WVDWSB',wvdwsb,0.0d0)
801 call reada(weightcard,'WELSB',welsb,0.0D0)
802 call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
803 call reada(weightcard,'WANG_NUCL',wang_nucl,0.0D0)
804 call reada(weightcard,'WSBLOC',wsbloc,0.0D0)
805 call reada(weightcard,'WTOR_NUCL',wtor_nucl,0.0D0)
806 ! print *,"KUR..",wtor_nucl
807 call reada(weightcard,'WTORD_NUCL',wtor_d_nucl,0.0D0)
808 call reada(weightcard,'WCORR_NUCL',wcorr_nucl,0.0D0)
809 call reada(weightcard,'WCORR3_NUC',wcorr3_nucl,0.0D0)
810 call reada(weightcard,'WBOND',wbond,1.0D0)
811 call reada(weightcard,'WTOR',wtor,1.0D0)
812 call reada(weightcard,'WTORD',wtor_d,1.0D0)
813 call reada(weightcard,'WSHIELD',wshield,0.05D0)
814 call reada(weightcard,'LIPSCALE',lipscale,1.0D0)
815 call reada(weightcard,'WLT',wliptran,1.0D0)
816 call reada(weightcard,'WTUBE',wtube,1.0d0)
817 call reada(weightcard,'WANG',wang,1.0D0)
818 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
819 call reada(weightcard,'SCAL14',scal14,0.4D0)
820 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
821 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
822 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
823 call reada(weightcard,'TEMP0',temp0,300.0d0)
824 call reada(weightcard,'WSCBASE',wscbase,0.0D0)
825 if (index(weightcard,'SOFT').gt.0) ipot=6
826 call reada(weightcard,'WBOND_NUCL',wbond_nucl,0.0D0)
827 call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
828 call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
829 call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
830 call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
831 call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
832 call reada(weightcard,'WSCPHO',wscpho,0.0d0)
833 call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
835 ! 12/1/95 Added weight for the multi-body term WCORR
836 call reada(weightcard,'WCORRH',wcorr,1.0D0)
837 if (wcorr4.gt.0.0d0) wcorr=wcorr4
857 weights(26)=wvdwpp_nucl
863 weights(32)=wbond_nucl
864 weights(33)=wang_nucl
866 weights(35)=wtor_nucl
867 weights(36)=wtor_d_nucl
868 weights(37)=wcorr_nucl
869 weights(38)=wcorr3_nucl
877 weights(56)=wcat_tran
878 if(me.eq.king.or..not.out1file) &
879 write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
880 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
882 10 format (/'Energy-term weights (unscaled):'// &
883 'WSCC= ',f10.6,' (SC-SC)'/ &
884 'WSCP= ',f10.6,' (SC-p)'/ &
885 'WELEC= ',f10.6,' (p-p electr)'/ &
886 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
887 'WBOND= ',f10.6,' (stretching)'/ &
888 'WANG= ',f10.6,' (bending)'/ &
889 'WSCLOC= ',f10.6,' (SC local)'/ &
890 'WTOR= ',f10.6,' (torsional)'/ &
891 'WTORD= ',f10.6,' (double torsional)'/ &
892 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
893 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
894 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
895 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
896 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
897 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ &
898 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
899 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
900 'WTURN6= ',f10.6,' (turns, 6th order)')
901 if(me.eq.king.or..not.out1file)then
902 if (wcorr4.gt.0.0d0) then
903 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',&
904 'between contact pairs of peptide groups'
905 write (iout,'(2(a,f5.3/))') &
906 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,&
907 'Range of quenching the correlation terms:',2*delt_corr
908 else if (wcorr.gt.0.0d0) then
909 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',&
910 'between contact pairs of peptide groups'
912 write (iout,'(a,f8.3)') &
913 'Scaling factor of 1,4 SC-p interactions:',scal14
914 write (iout,'(a,f8.3)') &
915 'General scaling factor of SC-p interactions:',scalscp
917 r0_corr=cutoff_corr-delt_corr
919 aad(i,1)=scalscp*aad(i,1)
920 aad(i,2)=scalscp*aad(i,2)
921 bad(i,1)=scalscp*bad(i,1)
922 bad(i,2)=scalscp*bad(i,2)
924 call rescale_weights(t_bath)
925 if(me.eq.king.or..not.out1file) &
926 write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
927 wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
929 22 format (/'Energy-term weights (scaled):'// &
930 'WSCC= ',f10.6,' (SC-SC)'/ &
931 'WSCP= ',f10.6,' (SC-p)'/ &
932 'WELEC= ',f10.6,' (p-p electr)'/ &
933 'WVDWPP= ',f10.6,' (p-p VDW)'/ &
934 'WBOND= ',f10.6,' (stretching)'/ &
935 'WANG= ',f10.6,' (bending)'/ &
936 'WSCLOC= ',f10.6,' (SC local)'/ &
937 'WTOR= ',f10.6,' (torsional)'/ &
938 'WTORD= ',f10.6,' (double torsional)'/ &
939 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/ &
940 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/ &
941 'WCORR4= ',f10.6,' (multi-body 4th order)'/ &
942 'WCORR5= ',f10.6,' (multi-body 5th order)'/ &
943 'WCORR6= ',f10.6,' (multi-body 6th order)'/ &
944 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ &
945 'WTURN3= ',f10.6,' (turns, 3rd order)'/ &
946 'WTURN4= ',f10.6,' (turns, 4th order)'/ &
947 'WTURN6= ',f10.6,' (turns, 6th order)')
948 if(me.eq.king.or..not.out1file) &
949 write (iout,*) "Reference temperature for weights calculation:",&
951 call reada(weightcard,"D0CM",d0cm,3.78d0)
952 call reada(weightcard,"AKCM",akcm,15.1d0)
953 call reada(weightcard,"AKTH",akth,11.0d0)
954 call reada(weightcard,"AKCT",akct,12.0d0)
955 call reada(weightcard,"V1SS",v1ss,-1.08d0)
956 call reada(weightcard,"V2SS",v2ss,7.61d0)
957 call reada(weightcard,"V3SS",v3ss,13.7d0)
958 call reada(weightcard,"EBR",ebr,-5.50D0)
959 call reada(weightcard,"ATRISS",atriss,0.301D0)
960 call reada(weightcard,"BTRISS",btriss,0.021D0)
961 call reada(weightcard,"CTRISS",ctriss,1.001D0)
962 call reada(weightcard,"DTRISS",dtriss,1.001D0)
963 call reada(weightcard,"SSSCALE",ssscale,1.0D0)
964 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
966 call reada(weightcard,"HT",Ht,0.0D0)
968 ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
969 Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
970 akcm=akcm/wsc*ssscale
971 akth=akth/wsc*ssscale
972 akct=akct/wsc*ssscale
973 v1ss=v1ss/wsc*ssscale
974 v2ss=v2ss/wsc*ssscale
975 v3ss=v3ss/wsc*ssscale
977 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
980 if(me.eq.king.or..not.out1file) then
981 write (iout,*) "Parameters of the SS-bond potential:"
982 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,&
984 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
985 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
986 write (iout,*)" HT",Ht
987 print *,'indpdb=',indpdb,' pdbref=',pdbref
989 if (indpdb.gt.0 .or. pdbref) then
990 read(inp,'(a)') pdbfile
991 if(me.eq.king.or..not.out1file) &
992 write (iout,'(2a)') 'PDB data will be read from file ',&
993 pdbfile(:ilen(pdbfile))
994 open(ipdbin,file=pdbfile,status='old',err=33)
996 33 write (iout,'(a)') 'Error opening PDB file.'
999 ! print *,'Begin reading pdb data'
1001 if (.not.allocated(crefjlee)) allocate (crefjlee(3,2*nres+2))
1004 crefjlee(j,i)=c(j,i)
1009 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
1010 & (crefjlee(j,i+nres),j=1,3)
1014 ! call int_from_cart1(.true.)
1016 ! print *,'Finished reading pdb data'
1017 if(me.eq.king.or..not.out1file) &
1018 write (iout,'(a,i7,a,i7)')'nsup=',nsup,&
1019 ' nstart_sup=',nstart_sup !,"ergwergewrgae"
1020 !el if(.not.allocated(itype_pdb))
1021 allocate(itype_pdb(nres))
1023 itype_pdb(i)=itype(i,1)
1027 nct=nstart_sup+nsup-1
1028 !el if(.not.allocated(icont_ref))
1029 allocate(icont_ref(2,(nres/2)*nres)) ! maxcont=12*maxres
1030 call contact(.false.,ncont_ref,icont_ref,co)
1033 if(me.eq.king.or..not.out1file) &
1034 write(iout,*)'Adding sidechains'
1038 if (iti.ne.10 .and. itype(i,1).ne.ntyp1) then
1041 do while (fail.and.nsi.le.maxsi)
1042 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail,molnum(i))
1045 if(fail) write(iout,*)'Adding sidechain failed for res ',&
1046 i,' after ',nsi,' trials'
1051 print *,"CZY TU DOCHODZE"
1052 if (indpdb.eq.0) then
1054 allocate(sequence(maxres,5))
1058 ! Read sequence if not taken from the pdb file.
1060 read (inp,*) nres_molec(molec)
1061 print *,'nres=',nres
1062 if (iscode.gt.0) then
1063 read (inp,'(80a1)') (sequence(i,molec)(1:1),i=1,nres_molec(molec))
1065 read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1067 ! read(inp,*) weightcard_t
1068 ! print *,"po seq" weightcard_t
1069 ! Convert sequence to numeric code
1071 do i=1,nres_molec(molec)
1073 itype(i,1)=rescode(i,sequence(i,molec),iscode,molec)
1078 ! read(inp,*) weightcard_t
1079 ! print *,"po seq", weightcard_t
1082 ! Read sequence if not taken from the pdb file.
1084 read (inp,*) nres_molec(molec)
1085 print *,'nres=',nres_molec(molec)
1086 ! allocate(sequence(maxres,5))
1087 ! if (iscode.gt.0) then
1088 read (inp,'(20a4)') (sequence(i,molec),i=1,nres_molec(molec))
1090 print *,(sequence(i,molec),i=1,nres_molec(molec))
1091 ! Convert sequence to numeric code
1093 do i=1,nres_molec(molec)
1095 istype(itmp)=sugarcode(sequence(i,molec)(1:1),i)
1096 sequence(i,molec)=sequence(i,molec)(1:2)
1097 itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1098 write(iout,*) "NUCLE=", itype(itmp,molec)
1103 ! Read sequence if not taken from the pdb file.
1105 read (inp,*) nres_molec(molec)
1106 ! print *,'nres=',nres
1107 read (inp,'(20(1x,a3))') (sequence(i,molec),i=1,nres_molec(molec))
1108 ! Convert sequence to numeric code
1109 print *,nres_molec(molec)
1110 do i=1,nres_molec(molec)
1113 itype(itmp,molec)=rescode(i,sequence(i,molec),iscode,molec)
1118 nres=nres+nres_molec(i)
1119 print *,"nres_molec",nres,nres_molec(i)
1122 ! Assign initial virtual bond lengths
1123 if(.not.allocated(molnum)) then
1124 allocate(molnum(nres+1))
1127 do j=1,nres_molec(i)
1132 ! print *,nres_molec(i)
1135 if(.not.allocated(vbld)) then
1136 print *, "I DO ENTER"
1137 allocate(vbld(2*nres))
1139 if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
1142 if (molnum(i).eq.1) then
1147 write(iout,*) "typy",itype(i,mnum),ntyp1_molec(mnum),i
1150 if ((itype(i,mnum).eq.ntyp1_molec(mnum)).or.&
1151 (itype(i-1,mnum).eq.ntyp1_molec(mnum))) then
1159 if (molnum(i).eq.1) then
1160 ! print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
1161 vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
1162 vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
1164 write(iout,*) "typy2",itype(i,mnum),ntyp1_molec(mnum),i
1165 if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
1166 vbld_inv(i+nres)=1.0
1169 vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1170 vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
1173 ! write (iout,*) "i",i," itype",itype(i,1),
1174 ! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
1178 ! print '(20i4)',(itype(i,1),i=1,nres)
1179 !----------------------------
1180 !el reallocate tables
1183 ! c_alloc(j,i)=c(j,i)
1184 ! dc_alloc(j,i)=dc(j,i)
1188 !elwrite(iout,*) "itype",i,itype(i,1)
1189 ! itype_alloc(i)=itype(i,1)
1195 ! allocate(c(3,2*nres+4))
1196 ! allocate(dc(3,0:2*nres+2))
1197 ! allocate(itype(nres+2))
1198 allocate(itel(nres+2))
1203 ! c(j,i)=c_alloc(j,i)
1204 ! dc(j,i)=dc_alloc(j,i)
1208 ! itype(i,1)=itype_alloc(i)
1211 !--------------------------
1214 if (itype(i,1).eq.ntyp1 .or. itype(i+1,1).eq.ntyp1) then
1216 if (itype(i,1).eq.ntyp1) then
1220 else if (iabs(itype(i+1,1)).ne.20) then
1222 else if (iabs(itype(i,1)).ne.20) then
1229 if(me.eq.king.or..not.out1file)then
1230 write (iout,*) "ITEL"
1233 write (iout,*) i,itype(i,1),itel(i)
1235 print *,'Call Read_Bridge.'
1238 !--------------------------------
1239 ! print *,"tu dochodze"
1240 ! znamy nres oraz nss można zaalokowac potrzebne tablice
1241 call alloc_geo_arrays
1242 print *,"after alloc_geo_arrays"
1243 call alloc_ener_arrays
1244 print *,"after alloc_ener_arrays"
1245 !--------------------------------
1246 ! 8/13/98 Set limits to generating the dihedral angles
1251 read (inp,*) ndih_constr
1252 if (ndih_constr.gt.0) then
1254 allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr)
1255 allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr)
1256 allocate(ftors(ndih_constr)) !(maxdih_constr)
1258 ! read (inp,*) ftors
1259 read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), &
1261 if(me.eq.king.or..not.out1file)then
1263 'There are',ndih_constr,' constraints on phi angles.'
1265 write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), &
1270 phi0(i)=deg2rad*phi0(i)
1271 drange(i)=deg2rad*drange(i)
1273 ! if(me.eq.king.or..not.out1file) &
1274 ! write (iout,*) 'FTORS',ftors
1277 phibound(1,ii) = phi0(i)-drange(i)
1278 phibound(2,ii) = phi0(i)+drange(i)
1280 else if (ndih_constr.lt.0) then
1282 allocate(idih_constr(nres))
1283 allocate(secprob(3,nres))
1284 allocate(vpsipred(3,nres))
1285 allocate(sdihed(2,nres))
1286 call card_concat(weightcard,.true.)
1287 call reada(weightcard,"PHIHEL",phihel,50.0D0)
1288 call reada(weightcard,"PHIBET",phibet,180.0D0)
1289 call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0)
1290 call reada(weightcard,"SIGMABET",sigmabet,40.0d0)
1291 call reada(weightcard,"WDIHC",wdihc,0.591D0)
1292 write (iout,*) "Weight of dihedral angle restraints",wdihc
1293 read(inp,'(9x,3f7.3)') &
1294 (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct)
1295 write (iout,*) "The secprob array"
1297 write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3)
1301 if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 &
1302 .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then
1303 ndih_constr=ndih_constr+1
1304 idih_constr(ndih_constr)=i
1307 vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2)
1308 sumv=sumv+vpsipred(j,ndih_constr)
1311 vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv
1313 phibound(1,ndih_constr)=phihel*deg2rad
1314 phibound(2,ndih_constr)=phibet*deg2rad
1315 sdihed(1,ndih_constr)=sigmahel*deg2rad
1316 sdihed(2,ndih_constr)=sigmabet*deg2rad
1321 if (with_theta_constr) then
1322 !C with_theta_constr is keyword allowing for occurance of theta constrains
1323 read (inp,*) ntheta_constr
1324 !C ntheta_constr is the number of theta constrains
1325 if (ntheta_constr.gt.0) then
1326 !C read (inp,*) ftors
1327 allocate(itheta_constr(ntheta_constr))
1328 allocate(theta_constr0(ntheta_constr))
1329 allocate(theta_drange(ntheta_constr),for_thet_constr(ntheta_constr))
1330 read (inp,*) (itheta_constr(i),theta_constr0(i), &
1331 theta_drange(i),for_thet_constr(i), &
1333 !C the above code reads from 1 to ntheta_constr
1334 !C itheta_constr(i) residue i for which is theta_constr
1335 !C theta_constr0 the global minimum value
1336 !C theta_drange is range for which there is no energy penalty
1337 !C for_thet_constr is the force constant for quartic energy penalty
1339 if(me.eq.king.or..not.out1file)then
1341 'There are',ntheta_constr,' constraints on phi angles.'
1342 do i=1,ntheta_constr
1343 write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), &
1348 do i=1,ntheta_constr
1349 theta_constr0(i)=deg2rad*theta_constr0(i)
1350 theta_drange(i)=deg2rad*theta_drange(i)
1352 !C if(me.eq.king.or..not.out1file)
1353 !C & write (iout,*) 'FTORS',ftors
1354 !C do i=1,ntheta_constr
1355 !C ii = itheta_constr(i)
1356 !C thetabound(1,ii) = phi0(i)-drange(i)
1357 !C thetabound(2,ii) = phi0(i)+drange(i)
1359 endif ! ntheta_constr.gt.0
1360 endif! with_theta_constr
1364 if (me.eq.king) then
1366 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1368 write (iout,'(a3,i5,2f10.1)') &
1369 restyp(itype(i,1),1),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1375 allocate(ireschain(nres))
1377 write(iout,*),"before seq2chains",ireschain
1379 write(iout,*) "after seq2chains",nchain
1380 allocate ( chain_border1(2,nchain))
1381 chain_border1(1,1)=1
1382 chain_border1(2,1)=chain_border(2,1)+1
1384 chain_border1(1,i)=chain_border(1,i)-1
1385 chain_border1(2,i)=chain_border(2,i)+1
1387 if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
1388 chain_border1(2,nchain)=nres
1389 write(iout,*) "nres",nres," nchain",nchain
1391 write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
1392 chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
1394 allocate(tabpermchain(50,5040))
1395 call chain_symmetry(npermchain,tabpermchain)
1396 print *,'NNT=',NNT,' NCT=',NCT
1397 if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
1398 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
1400 if(me.eq.king.or..not.out1file) &
1401 write (iout,'(a,i7)') 'nsup=',nsup
1403 if (nsup.le.(nct-nnt+1)) then
1404 do i=0,nct-nnt+1-nsup
1405 if (seq_comp(itype(nnt+i,1),itype_pdb(nstart_sup),nsup)) then
1410 write (iout,'(a)') &
1411 'Error - sequences to be superposed do not match.'
1414 do i=0,nsup-(nct-nnt+1)
1415 if (seq_comp(itype(nnt,1),itype_pdb(nstart_sup+i),nct-nnt+1)) &
1417 nstart_sup=nstart_sup+i
1422 write (iout,'(a)') &
1423 'Error - sequences to be superposed do not match.'
1426 if (nsup.eq.0) nsup=nct-nnt+1
1427 if (nstart_sup.eq.0) nstart_sup=nnt
1428 if (nstart_seq.eq.0) nstart_seq=nnt
1429 if(me.eq.king.or..not.out1file) &
1430 write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,&
1431 ' nstart_seq=',nstart_seq !,"242343453254"
1433 !--- Zscore rms -------
1434 if (nz_start.eq.0) nz_start=nnt
1435 if (nz_end.eq.0 .and. nsup.gt.0) then
1437 else if (nz_end.eq.0) then
1440 if(me.eq.king.or..not.out1file)then
1441 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1442 write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct
1444 !----------------------
1447 if (.not.pdbref) then
1448 call read_angles(inp,*38)
1450 38 write (iout,'(a)') 'Error reading reference structure.'
1452 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1453 stop 'Error reading reference structure'
1457 !zscore call geom_to_var(nvar,coord_exp_zs(1,1))
1464 cref(j,i,kkk)=c(j,i)
1467 call contact(.true.,ncont_ref,icont_ref,co)
1469 ! write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1471 !EL if (constr_dist.gt.0) call read_dist_constr
1472 !EL write (iout,*) "After read_dist_constr nhpb",nhpb
1473 !EL if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1474 !EL call hpb_partition
1475 if(me.eq.king.or..not.out1file) &
1476 write (iout,*) 'Contact order:',co
1478 if(me.eq.king.or..not.out1file) &
1479 write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1482 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1484 if(me.eq.king.or..not.out1file) &
1485 write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1486 icont_ref(1,i),' ',&
1487 restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1490 if (constr_homology.gt.0) then
1491 ! write (iout,*) "Calling read_constr_homology"
1493 call read_constr_homology
1494 if (indpdb.gt.0 .or. pdbref) then
1497 c(j,i)=crefjlee(j,i)
1498 cref(j,i,1)=crefjlee(j,i)
1504 write (iout,*) "sc_loc_geom: Array C"
1506 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
1509 write (iout,*) "Array Cref"
1511 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
1512 (cref(j,i+nres,1),j=1,3)
1515 call int_from_cart1(.false.)
1516 call sc_loc_geom(.false.)
1518 thetaref(i)=theta(i)
1523 dc(j,i)=c(j,i+1)-c(j,i)
1524 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1529 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1530 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1535 if (start_from_model) then
1538 read(inp,'(a)',end=332,err=332) pdbfile
1539 if (me.eq.king .or. .not. out1file)&
1540 write (iout,'(a,5x,a)') 'Opening PDB file',&
1541 pdbfile(:ilen(pdbfile))
1542 open(ipdbin,file=pdbfile,status='old',err=336)
1544 336 write (iout,'(a,5x,a)') 'Error opening PDB file',&
1545 pdbfile(:ilen(pdbfile))
1552 call readpdb_template(nmodel_start+1)
1554 if (nres.ge.nres_temp) then
1555 nmodel_start=nmodel_start+1
1556 pdbfiles_chomo(nmodel_start)=pdbfile
1559 chomo(j,i,nmodel_start)=c(j,i)
1563 if (me.eq.king .or. .not. out1file) &
1564 write (iout,'(a,2i7,1x,a)') &
1565 "Different number of residues",nres_temp,nres, &
1571 if (nmodel_start.eq.0) then
1572 if (me.eq.king .or. .not. out1file) &
1573 write (iout,'(a)') &
1574 "No valid starting model found START_FROM_MODELS is OFF"
1575 start_from_model=.false.
1577 write (iout,*) "nmodel_start",nmodel_start
1582 if (constr_dist.gt.0) call read_dist_constr
1583 write (iout,*) "After read_dist_constr nhpb",nhpb
1584 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1585 if (velnanoconst.ne.0) call read_afmnano
1588 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1589 .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1590 modecalc.ne.10) then
1591 ! If input structure hasn't been supplied from the PDB file read or generate
1593 if (iranconf.eq.0 .and. .not. extconf) then
1594 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1595 write (iout,'(a)') 'Initial geometry will be read in.'
1597 read(inp,'(8f10.5)',end=36,err=36) &
1598 ((c(l,k),l=1,3),k=1,nres),&
1599 ((c(l,k+nres),l=1,3),k=nnt,nct)
1600 write (iout,*) "Exit READ_CART"
1601 write (iout,'(8f10.5)') &
1602 ((c(l,k),l=1,3),k=1,nres)
1603 write (iout,'(8f10.5)') &
1604 ((c(l,k+nres),l=1,3),k=nnt,nct)
1605 call int_from_cart1(.true.)
1606 write (iout,*) "Finish INT_TO_CART"
1609 dc(j,i)=c(j,i+1)-c(j,i)
1610 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1614 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1616 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1617 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1623 write(iout,*) "read angles from input"
1624 call read_angles(inp,*36)
1629 36 write (iout,'(a)') 'Error reading angle file.'
1631 call mpi_finalize( MPI_COMM_WORLD,IERR )
1633 stop 'Error reading angle file.'
1635 else if (extconf) then
1636 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1637 write (iout,'(a)') 'Extended chain initial geometry.'
1639 theta(i)=90d0*deg2rad
1640 if (molnum(i).eq.2) theta(i)=160d0*deg2rad
1643 phi(i)=180d0*deg2rad
1644 if (molnum(i).eq.2) phi(i)=30d0*deg2rad
1647 alph(i)=110d0*deg2rad
1648 if (molnum(i).eq.2) alph(i)=30d0*deg2rad
1651 omeg(i)=-120d0*deg2rad
1652 if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
1653 if (itype(i,1).le.0) omeg(i)=-omeg(i)
1657 if(me.eq.king.or..not.out1file) &
1658 write (iout,'(a)') 'Random-generated initial geometry.'
1662 if (me.eq.king .or. fg_rank.eq.0 .and. &
1663 ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1667 call gen_rand_conf(itmp,*30)
1669 30 write (iout,*) 'Failed to generate random conformation',&
1671 write (*,*) 'Processor:',me,&
1672 ' Failed to generate random conformation',&
1682 write (iout,'(a,i3,a)') 'Processor:',me,&
1683 ' error in generating random conformation.'
1684 write (*,'(a,i3,a)') 'Processor:',me,&
1685 ' error in generating random conformation.'
1688 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1694 call gen_rand_conf(itmp,*335)
1696 335 write (iout,*) 'Failed to generate random conformation',&
1698 write (*,*) 'Failed to generate random conformation',&
1701 write (iout,'(a,i3,a)') 'Processor:',me,&
1702 ' error in generating random conformation.'
1703 write (*,'(a,i3,a)') 'Processor:',me,&
1704 ' error in generating random conformation.'
1709 elseif (modecalc.eq.4) then
1710 read (inp,'(a)') intinname
1711 open (intin,file=intinname,status='old',err=333)
1712 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1713 write (iout,'(a)') 'intinname',intinname
1714 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1716 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1718 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1720 stop 'Error opening angle file.'
1724 ! Generate distance constraints, if the PDB structure is to be regularized.
1725 if (nthread.gt.0) then
1726 call read_threadbase
1729 if (me.eq.king .or. .not. out1file) &
1731 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1732 write (iout,'(/a,i3,a)') &
1733 'The chain contains',ns,' disulfide-bridging cysteines.'
1734 write (iout,'(20i4)') (iss(i),i=1,ns)
1736 write(iout,*)"Running with dynamic disulfide-bond formation"
1738 write (iout,'(/a/)') 'Pre-formed links are:'
1744 if (me.eq.king.or..not.out1file) &
1745 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1746 restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1752 if (ns.gt.0.and.dyn_ss) then
1756 forcon(i-nss)=forcon(i)
1763 dyn_ss_mask(iss(i))=.true.
1766 if (i2ndstr.gt.0) call secstrp2dihc
1767 if (indpdb.gt.0) then
1768 write(iout,*) "WCHODZE TU!!"
1769 call int_from_cart1(.true.)
1771 ! call geom_to_var(nvar,x)
1772 ! call etotal(energia(0))
1773 ! call enerprint(energia(0))
1774 ! call briefout(0,etot)
1776 !d write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1777 !d write (iout,'(a)') 'Variable list:'
1778 !d write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1780 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1781 write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1782 'Processor',myrank,': end reading molecular data.'
1785 end subroutine molread
1786 !-----------------------------------------------------------------------------
1787 subroutine read_constr_homology
1788 use control, only:init_int_table,homology_partition
1789 use MD_data, only:iset
1791 ! include 'DIMENSIONS'
1795 ! include 'COMMON.SETUP'
1796 ! include 'COMMON.CONTROL'
1797 ! include 'COMMON.HOMOLOGY'
1798 ! include 'COMMON.CHAIN'
1799 ! include 'COMMON.IOUNITS'
1800 ! include 'COMMON.MD'
1801 ! include 'COMMON.QRESTR'
1802 ! include 'COMMON.GEO'
1803 ! include 'COMMON.INTERACT'
1804 ! include 'COMMON.NAMES'
1805 ! include 'COMMON.VAR'
1808 ! double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1810 ! common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1811 ! & sigma_odl_temp(maxres,maxres,max_template)
1813 character*24 model_ki_dist, model_ki_angle
1814 character*500 controlcard
1815 integer :: ki,i,ii,j,k,l
1816 integer, dimension (:), allocatable :: ii_in_use
1817 integer :: i_tmp,idomain_tmp,&
1818 irec,ik,iistart,nres_temp
1821 logical :: liiflag,lfirst
1824 ! FP - Nov. 2014 Temporary specifications for new vars
1826 real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1827 rescore3_tmp, dist_cut
1828 real(kind=8), dimension (:,:),allocatable :: rescore
1829 real(kind=8), dimension (:,:),allocatable :: rescore2
1830 real(kind=8), dimension (:,:),allocatable :: rescore3
1831 real(kind=8) :: distal
1832 character*24 tpl_k_rescore
1833 character*256 pdbfile
1835 ! -----------------------------------------------------------------
1836 ! Reading multiple PDB ref structures and calculation of retraints
1837 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1839 ! -----------------------------------------------------------------
1842 ! Alternative: reading from input
1843 call card_concat(controlcard,.true.)
1844 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1845 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1846 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1847 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1848 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1849 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1850 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1851 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1852 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1853 if(.not.read2sigma.and.start_from_model) then
1854 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1855 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1856 start_from_model=.false.
1857 iranconf=(indpdb.le.0)
1859 if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1860 write(iout,*) 'START_FROM_MODELS is ON'
1861 ! if(start_from_model .and. rest) then
1862 ! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1863 ! write(iout,*) 'START_FROM_MODELS is OFF'
1864 ! write(iout,*) 'remove restart keyword from input'
1867 if (start_from_model) nmodel_start=constr_homology
1868 if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1869 if (homol_nset.gt.1)then
1870 call card_concat(controlcard,.true.)
1871 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1872 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1873 ! write(iout,*) "iset homology_weight "
1875 write(iout,*) i,waga_homology(i)
1878 iset=mod(kolor,homol_nset)+1
1881 waga_homology(1)=1.0
1884 !d write (iout,*) "nnt",nnt," nct",nct
1887 if (read_homol_frag) then
1888 call read_klapaucjusz
1894 ! write(iout,*) 'nnt=',nnt,'nct=',nct
1897 ! do k=1,constr_homology
1911 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1912 if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1914 do k=1,constr_homology
1916 read(inp,'(a)') pdbfile
1917 pdbfiles_chomo(k)=pdbfile
1918 if(me.eq.king .or. .not. out1file) &
1919 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1920 pdbfile(:ilen(pdbfile))
1921 open(ipdbin,file=pdbfile,status='old',err=33)
1923 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
1924 pdbfile(:ilen(pdbfile))
1927 ! print *,'Begin reading pdb data'
1929 ! Files containing res sim or local scores (former containing sigmas)
1932 write(kic2,'(bz,i2.2)') k
1934 tpl_k_rescore="template"//kic2//".sco"
1935 write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1938 write(iout,*) "read2sigma",read2sigma
1940 if (read2sigma) then
1941 call readpdb_template(k)
1945 write(iout,*) "after readpdb"
1946 if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1949 if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1950 if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1951 if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1952 if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1953 if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1954 if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1955 if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1956 if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1957 if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1958 if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1959 if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1960 if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1961 if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1962 if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1963 ! if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1964 if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1965 if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1966 if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1967 if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1968 ! if(.not.allocated(distance)) allocate(distance(constr_homology))
1969 ! if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1973 ! Distance restraints
1976 ! Copy the coordinates from reference coordinates (?)
1977 do i=1,2*nres_chomo(k)
1980 ! write (iout,*) "c(",j,i,") =",c(j,i)
1984 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1986 ! write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1987 open (ientin,file=tpl_k_rescore,status='old')
1988 if (nnt.gt.1) rescore(k,1)=0.0d0
1989 do irec=nnt,nct ! loop for reading res sim
1990 if (read2sigma) then
1991 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1992 rescore3_tmp,idomain_tmp
1994 idomain(k,i_tmp)=idomain_tmp
1995 rescore(k,i_tmp)=rescore_tmp
1996 rescore2(k,i_tmp)=rescore2_tmp
1997 rescore3(k,i_tmp)=rescore3_tmp
1998 if (.not. out1file .or. me.eq.king)&
1999 write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
2000 i_tmp,rescore2_tmp,rescore_tmp,&
2001 rescore3_tmp,idomain_tmp
2004 read (ientin,*,end=1401) rescore_tmp
2006 ! rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2007 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2008 ! write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2013 if (waga_dist.ne.0.0d0) then
2021 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2022 ! write (iout,*) k,i,j,distal,dist2_cut
2024 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2025 .and. distal.le.dist2_cut ) then
2031 ! write (iout,*) "k",k
2032 ! write (iout,*) "i",i," j",j," constr_homology",
2037 if (read2sigma) then
2040 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2042 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2043 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2044 sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2046 if (odl(k,ii).le.dist_cut) then
2047 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2050 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2051 dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2053 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2054 dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2058 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2061 ! l_homo(k,ii)=.false.
2067 ! write (iout,*) "Distance restraints set"
2070 ! Theta, dihedral and SC retraints
2072 if (waga_angle.gt.0.0d0) then
2073 ! open (ientin,file=tpl_k_sigma_dih,status='old')
2074 ! do irec=1,maxres-3 ! loop for reading sigma_dih
2075 ! read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2076 ! if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2077 ! sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2078 ! & sigma_dih(k,i+nnt-1)
2083 if (idomain(k,i).eq.0) then
2087 dih(k,i)=phiref(i) ! right?
2088 ! read (ientin,*) sigma_dih(k,i) ! original variant
2089 ! write (iout,*) "dih(",k,i,") =",dih(k,i)
2090 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2091 ! & "rescore(",k,i-1,") =",rescore(k,i-1),
2092 ! & "rescore(",k,i-2,") =",rescore(k,i-2),
2093 ! & "rescore(",k,i-3,") =",rescore(k,i-3)
2095 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2096 rescore(k,i-2)+rescore(k,i-3))/4.0
2097 ! if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2098 ! write (iout,*) "Raw sigmas for dihedral angle restraints"
2099 ! write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2100 ! sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2101 ! rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2102 ! Instead of res sim other local measure of b/b str reliability possible
2103 if (sigma_dih(k,i).ne.0) &
2104 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2105 ! sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2109 ! write (iout,*) "Dihedral angle restraints set"
2112 if (waga_theta.gt.0.0d0) then
2113 ! open (ientin,file=tpl_k_sigma_theta,status='old')
2114 ! do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2115 ! read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2116 ! sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2117 ! & sigma_theta(k,i+nnt-1)
2122 do i = nnt+2,nct ! right? without parallel.
2123 ! do i = i=1,nres ! alternative for bounds acc to readpdb?
2124 ! do i=ithet_start,ithet_end ! with FG parallel.
2125 if (idomain(k,i).eq.0) then
2126 sigma_theta(k,i)=0.0
2129 thetatpl(k,i)=thetaref(i)
2130 ! write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2131 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2132 ! & "rescore(",k,i-1,") =",rescore(k,i-1),
2133 ! & "rescore(",k,i-2,") =",rescore(k,i-2)
2134 ! read (ientin,*) sigma_theta(k,i) ! 1st variant
2135 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2137 ! if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2138 if (sigma_theta(k,i).ne.0) &
2139 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2141 ! sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2142 ! rescore(k,i-2) ! right expression ?
2143 ! sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2146 ! write (iout,*) "Angle restraints set"
2149 if (waga_d.gt.0.0d0) then
2150 ! open (ientin,file=tpl_k_sigma_d,status='old')
2151 ! do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2152 ! read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2153 ! sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2154 ! & sigma_d(k,i+nnt-1)
2158 do i = nnt,nct ! right? without parallel.
2159 ! do i=2,nres-1 ! alternative for bounds acc to readpdb?
2160 ! do i=loc_start,loc_end ! with FG parallel.
2161 if (itype(i,1).eq.10) cycle
2162 if (idomain(k,i).eq.0 ) then
2169 ! write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2170 ! write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2171 ! write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2172 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2173 sigma_d(k,i)=rescore3(k,i) ! right expression ?
2174 if (sigma_d(k,i).ne.0) &
2175 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2177 ! sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2178 ! sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2179 ! read (ientin,*) sigma_d(k,i) ! 1st variant
2183 ! write (iout,*) "SC restraints set"
2186 ! remove distance restraints not used in any model from the list
2187 ! shift data in all arrays
2189 ! write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
2190 if (waga_dist.ne.0.0d0) then
2197 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2198 ! & .and. distal.le.dist2_cut ) then
2199 ! write (iout,*) "i",i," j",j," ii",ii
2201 if (ii_in_use(ii).eq.0.and.liiflag.or. &
2202 ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2209 if(i10.eq.lim_odl) i10=i10+1
2211 ires_homo(iistart+ki)=ires_homo(ki+i01)
2212 jres_homo(iistart+ki)=jres_homo(ki+i01)
2213 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2214 do k=1,constr_homology
2215 odl(k,iistart+ki)=odl(k,ki+i01)
2216 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2217 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2220 iistart=iistart+i10-i01
2223 if (ii_in_use(ii).ne.0.and..not.liiflag) then
2231 ! write (iout,*) "Removing distances completed"
2233 endif ! .not. klapaucjusz
2235 if (constr_homology.gt.0) call homology_partition
2236 write (iout,*) "After homology_partition"
2238 if (constr_homology.gt.0) call init_int_table
2239 write (iout,*) "After init_int_table"
2241 ! endif ! .not. klapaucjusz
2243 ! if (constr_homology.gt.0) call homology_partition
2244 ! write (iout,*) "After homology_partition"
2246 ! if (constr_homology.gt.0) call init_int_table
2247 ! write (iout,*) "After init_int_table"
2249 ! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2250 ! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2254 if (.not.out_template_restr) return
2255 !d write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2256 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2257 write (iout,*) "Distance restraints from templates"
2259 write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
2260 ii,ires_homo(ii),jres_homo(ii),&
2261 (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
2262 ki=1,constr_homology)
2264 write (iout,*) "Dihedral angle restraints from templates"
2266 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2267 (rad2deg*dih(ki,i),&
2268 rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2270 write (iout,*) "Virtual-bond angle restraints from templates"
2272 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2273 (rad2deg*thetatpl(ki,i),&
2274 rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2276 write (iout,*) "SC restraints from templates"
2278 write(iout,'(i7,100(4f8.2,4x))') i,&
2279 (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
2280 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2284 end subroutine read_constr_homology
2285 !-----------------------------------------------------------------------------
2286 subroutine read_klapaucjusz
2289 ! include 'DIMENSIONS'
2293 ! include 'COMMON.SETUP'
2294 ! include 'COMMON.CONTROL'
2295 ! include 'COMMON.HOMOLOGY'
2296 ! include 'COMMON.CHAIN'
2297 ! include 'COMMON.IOUNITS'
2298 ! include 'COMMON.MD'
2299 ! include 'COMMON.GEO'
2300 ! include 'COMMON.INTERACT'
2301 ! include 'COMMON.NAMES'
2302 character(len=256) fragfile
2303 integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2305 integer, dimension(:,:), allocatable :: iresclust,inclust
2308 character(len=2) :: kic2
2309 character(len=24) :: model_ki_dist, model_ki_angle
2310 character(len=500) :: controlcard
2311 integer :: ki, i, j, jj,k, l, i_tmp,&
2313 ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2314 i01,i10,nnt_chain,nct_chain
2315 real(kind=8) :: distal
2316 logical :: lprn = .true.
2317 integer :: nres_temp
2320 logical :: liiflag,lfirst
2322 real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2323 real(kind=8), dimension (:,:), allocatable :: rescore
2324 real(kind=8), dimension (:,:), allocatable :: rescore2
2325 character(len=24) :: tpl_k_rescore
2326 character(len=256) :: pdbfile
2329 ! For new homol impl
2331 ! include 'COMMON.VAR'
2333 ! write (iout,*) "READ_KLAPAUCJUSZ"
2334 ! print *,"READ_KLAPAUCJUSZ"
2336 call getenv("FRAGFILE",fragfile)
2337 write (iout,*) "Opening", fragfile
2339 open(ientin,file=fragfile,status="old",err=10)
2340 ! write (iout,*) " opened"
2349 itype_temp(:)=itype(:,1)
2353 ! write (iout,*) "Entering loop"
2356 DO IGR = 1,NCHAIN_GROUP
2358 ! write (iout,*) "igr",igr
2360 read(ientin,*) constr_homology,nclust
2361 if (start_from_model) then
2362 nmodel_start=constr_homology
2369 ichain=iequiv(1,igr)
2370 nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2371 nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2372 ! write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2375 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2376 if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2377 do k=1,constr_homology
2378 read(ientin,'(a)') pdbfile
2379 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2380 pdbfile(:ilen(pdbfile))
2381 pdbfiles_chomo(k)=pdbfile
2382 open(ipdbin,file=pdbfile,status='old',err=33)
2384 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
2385 pdbfile(:ilen(pdbfile))
2390 call readpdb_template(k)
2400 read(ientin,*) ninclust(i),nresclust(i)
2401 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2402 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2405 ! Loop over clusters
2408 do ll = 1,ninclust(l)
2411 ! write (iout,*) "l",l," ll",ll," k",k
2417 idomain(k,iresclust(i,l)+1) = 1
2419 idomain(k,iresclust(i,l)) = 1
2423 ! Distance restraints
2426 ! Copy the coordinates from reference coordinates (?)
2432 ! write (iout,*) "c(",j,i,") =",c(j,i)
2435 call int_from_cart(.true.,.false.)
2436 call sc_loc_geom(.false.)
2438 thetaref(i)=theta(i)
2442 if (waga_dist.ne.0.0d0) then
2445 do i = nnt_chain,nct_chain-2
2452 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2453 ! write (iout,*) k,i,j,distal,dist2_cut
2455 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2456 .and. distal.le.dist2_cut ) then
2462 ! write (iout,*) "k",k
2463 ! write (iout,*) "i",i," j",j," constr_homology",
2465 ires_homo(ii)=i+chain_border1(1,igr)-1
2466 jres_homo(ii)=j+chain_border1(1,igr)-1
2468 if (read2sigma) then
2471 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2473 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2474 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2475 sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2477 if (odl(k,ii).le.dist_cut) then
2478 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2481 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2482 dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2484 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2485 dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2489 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2492 ! l_homo(k,ii)=.false.
2499 ! Theta, dihedral and SC retraints
2501 if (waga_angle.gt.0.0d0) then
2502 do i = nnt_chain+3,nct_chain
2503 iii=i+chain_border1(1,igr)-1
2504 if (idomain(k,i).eq.0) then
2505 ! sigma_dih(k,i)=0.0
2508 dih(k,iii)=phiref(i)
2510 (rescore(k,i)+rescore(k,i-1)+ &
2511 rescore(k,i-2)+rescore(k,i-3))/4.0
2512 ! write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2513 ! & " sigma_dihed",sigma_dih(k,i)
2514 if (sigma_dih(k,iii).ne.0) &
2515 sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2520 if (waga_theta.gt.0.0d0) then
2521 do i = nnt_chain+2,nct_chain
2522 iii=i+chain_border1(1,igr)-1
2523 if (idomain(k,i).eq.0) then
2524 ! sigma_theta(k,i)=0.0
2527 thetatpl(k,iii)=thetaref(i)
2528 sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2530 if (sigma_theta(k,iii).ne.0) &
2531 sigma_theta(k,iii)=1.0d0/ &
2532 (sigma_theta(k,iii)*sigma_theta(k,iii))
2536 if (waga_d.gt.0.0d0) then
2537 do i = nnt_chain,nct_chain
2538 iii=i+chain_border1(1,igr)-1
2539 if (itype(i,1).eq.10) cycle
2540 if (idomain(k,i).eq.0 ) then
2544 xxtpl(k,iii)=xxref(i)
2545 yytpl(k,iii)=yyref(i)
2546 zztpl(k,iii)=zzref(i)
2547 sigma_d(k,iii)=rescore(k,i)
2548 if (sigma_d(k,iii).ne.0) &
2549 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2550 ! if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2556 ! remove distance restraints not used in any model from the list
2557 ! shift data in all arrays
2559 ! write (iout,*) "ii_old",ii_old
2560 if (waga_dist.ne.0.0d0) then
2562 write (iout,*) "Distance restraints from templates"
2564 write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2565 iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2566 (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2567 ki=1,constr_homology)
2573 do i=nnt_chain,nct_chain-2
2576 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2577 ! & .and. distal.le.dist2_cut ) then
2578 ! write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
2580 if (ii_in_use(ii).eq.0.and.liiflag.or. &
2581 ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2588 if(i10.eq.lim_odl) i10=i10+1
2590 ires_homo(iistart+ki)=ires_homo(ki+i01)
2591 jres_homo(iistart+ki)=jres_homo(ki+i01)
2592 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2593 do k=1,constr_homology
2594 odl(k,iistart+ki)=odl(k,ki+i01)
2595 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2596 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2599 iistart=iistart+i10-i01
2602 if (ii_in_use(ii).ne.0.and..not.liiflag) then
2615 ichain=iequiv(i,igr)
2617 do j=nnt_chain,nct_chain
2618 jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
2619 do k=1,constr_homology
2621 sigma_dih(k,jj)=sigma_dih(k,j)
2622 thetatpl(k,jj)=thetatpl(k,j)
2623 sigma_theta(k,jj)=sigma_theta(k,j)
2624 xxtpl(k,jj)=xxtpl(k,j)
2625 yytpl(k,jj)=yytpl(k,j)
2626 zztpl(k,jj)=zztpl(k,j)
2627 sigma_d(k,jj)=sigma_d(k,j)
2631 jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
2632 ! write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
2633 do j=ii_old+1,lim_odl
2634 ires_homo(j+lll)=ires_homo(j)+jj
2635 jres_homo(j+lll)=jres_homo(j)+jj
2636 do k=1,constr_homology
2637 odl(k,j+lll)=odl(k,j)
2638 sigma_odl(k,j+lll)=sigma_odl(k,j)
2639 l_homo(k,j+lll)=l_homo(k,j)
2650 if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
2652 itype(:,1)=itype_temp(:)
2655 10 stop "Error in fragment file"
2656 end subroutine read_klapaucjusz
2658 !-----------------------------------------------------------
2659 subroutine seq2chains
2661 !c Split the total UNRES sequence, which has dummy residues separating
2662 !c the chains, into separate chains. The length of chain ichain is
2663 !c contained in chain_length(ichain), the first and last non-dummy
2664 !c residues are in chain_border(1,ichain) and chain_border(2,ichain),
2665 !c respectively. The lengths pertain to non-dummy residues only.
2668 ! include 'DIMENSIONS'
2669 use energy_data, only:molnum,nchain,chain_length,ireschain
2671 ! integer ireschain(nres)
2672 integer ii,ichain,i,j,mnum
2677 if (.not.allocated(chain_length)) allocate(chain_length(50))
2678 if (.not.allocated(chain_border)) allocate(chain_border(2,50))
2680 chain_length(ichain)=0
2682 do while (ii.lt.nres)
2683 write(iout,*) "in seq2chains",ii,new_chain
2685 if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
2686 if (.not.new_chain) then
2688 chain_border(2,ichain)=ii-1
2690 chain_border(1,ichain)=ii+1
2691 chain_length(ichain)=0
2695 chain_border(1,ichain)=ii
2698 chain_length(ichain)=chain_length(ichain)+1
2702 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
2705 chain_length(ichain)=chain_length(ichain)+1
2707 if (chain_length(ichain).gt.0) then
2708 chain_border(2,ichain)=ii
2715 do j=chain_border(1,i),chain_border(2,i)
2721 !---------------------------------------------------------------------
2722 subroutine chain_symmetry(npermchain,tabpermchain)
2724 !c Determine chain symmetry. nperm is the number of permutations and
2725 !c tabperchain contains the allowed permutations of the chains.
2728 ! include "DIMENSIONS"
2729 ! include "COMMON.IOUNITS"
2732 npermchain,tabpermchain(50,5040),&
2733 tabperm(50,5040),mapchain(50),&
2734 iequiv(50,nres),iflag(nres)
2735 integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
2736 nperm,npermc,ind,mnum
2737 if (nchain.eq.1) then
2740 !c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
2744 !c Look for equivalent chains
2746 write(iout,*) "nchain",nchain
2748 write(iout,*) "chain",i," from",chain_border(1,i),&
2749 " to",chain_border(2,i)
2751 "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
2759 if (iflag(i).gt.0) cycle
2761 nchain_group=nchain_group+1
2763 iequiv(iieq,nchain_group)=i
2765 if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
2767 !c do while(k.lt.chain_length(i) .and.
2768 !c & itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
2769 do k=0,chain_length(i)-1
2772 if (itype(chain_border(1,i)+k,mnum).ne.&
2773 itype(chain_border(1,j)+k,mnum)) exit
2775 if (k.lt.chain_length(i)) cycle
2778 iequiv(iieq,nchain_group)=j
2780 nequiv(nchain_group)=iieq
2782 write(iout,*) "Number of equivalent chain groups:",nchain_group
2783 write(iout,*) "Equivalent chain groups"
2785 write(iout,*) "group",i," #members",nequiv(i)," chains",&
2786 (iequiv(j,i),j=1,nequiv(i))
2792 mapchain(ind)=iequiv(j,i)
2795 write (iout,*) "mapchain"
2797 write (iout,*) i,mapchain(i)
2801 call permut(nequiv(i),nperm,tabperm)
2807 tabpermchain(k,j)=iequiv(tabperm(k,j),i)
2812 npermchain=npermchain*nperm
2818 tabpermchain(l,ind)=tabpermchain(l,j)
2821 tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
2830 itemp(mapchain(j))=tabpermchain(j,i)
2833 tabpermchain(j,i)=itemp(j)
2836 write(iout,*) "Number of chain permutations",npermchain
2837 write(iout,*) "Permutations"
2839 write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
2843 !c---------------------------------------------------------------------
2844 integer function tperm(i,iperm,tabpermchain)
2846 ! include 'DIMENSIONS'
2848 integer tabpermchain(50,5040)
2852 tperm=tabpermchain(i,iperm)
2857 !-----------------------------------------------------------------------------