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(.false.,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(.false.,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
1483 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1485 if(me.eq.king.or..not.out1file) &
1486 write (2,*) i,' ',restyp(itype(icont_ref(1,i),1),1),' ',&
1487 icont_ref(1,i),' ',&
1488 restyp(itype(icont_ref(2,i),1),1),' ',icont_ref(2,i)
1492 if (constr_homology.gt.0) then
1493 ! write (iout,*) "Calling read_constr_homology"
1495 call read_constr_homology
1496 if (indpdb.gt.0 .or. pdbref) then
1499 c(j,i)=crefjlee(j,i)
1500 cref(j,i,1)=crefjlee(j,i)
1506 write (iout,*) "sc_loc_geom: Array C"
1508 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),&
1511 write (iout,*) "Array Cref"
1513 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),&
1514 (cref(j,i+nres,1),j=1,3)
1517 call int_from_cart1(.false.)
1518 call sc_loc_geom(.false.)
1520 thetaref(i)=theta(i)
1525 dc(j,i)=c(j,i+1)-c(j,i)
1526 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1531 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1532 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1537 if (start_from_model) then
1540 read(inp,'(a)',end=332,err=332) pdbfile
1541 if (me.eq.king .or. .not. out1file)&
1542 write (iout,'(a,5x,a)') 'Opening PDB file',&
1543 pdbfile(:ilen(pdbfile))
1544 open(ipdbin,file=pdbfile,status='old',err=336)
1546 336 write (iout,'(a,5x,a)') 'Error opening PDB file',&
1547 pdbfile(:ilen(pdbfile))
1554 call readpdb_template(nmodel_start+1)
1556 if (nres.ge.nres_temp) then
1557 nmodel_start=nmodel_start+1
1558 pdbfiles_chomo(nmodel_start)=pdbfile
1561 chomo(j,i,nmodel_start)=c(j,i)
1565 if (me.eq.king .or. .not. out1file) &
1566 write (iout,'(a,2i7,1x,a)') &
1567 "Different number of residues",nres_temp,nres, &
1573 if (nmodel_start.eq.0) then
1574 if (me.eq.king .or. .not. out1file) &
1575 write (iout,'(a)') &
1576 "No valid starting model found START_FROM_MODELS is OFF"
1577 start_from_model=.false.
1579 write (iout,*) "nmodel_start",nmodel_start
1584 if (constr_dist.gt.0) call read_dist_constr
1585 write (iout,*) "After read_dist_constr nhpb",nhpb
1586 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
1587 if (velnanoconst.ne.0) call read_afmnano
1590 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
1591 .and. modecalc.ne.8 .and. modecalc.ne.9 .and. &
1592 modecalc.ne.10) then
1593 ! If input structure hasn't been supplied from the PDB file read or generate
1595 if (iranconf.eq.0 .and. .not. extconf) then
1596 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) &
1597 write (iout,'(a)') 'Initial geometry will be read in.'
1599 read(inp,'(8f10.5)',end=36,err=36) &
1600 ((c(l,k),l=1,3),k=1,nres),&
1601 ((c(l,k+nres),l=1,3),k=nnt,nct)
1602 write (iout,*) "Exit READ_CART"
1603 write (iout,'(8f10.5)') &
1604 ((c(l,k),l=1,3),k=1,nres)
1605 write (iout,'(8f10.5)') &
1606 ((c(l,k+nres),l=1,3),k=nnt,nct)
1607 call int_from_cart1(.true.)
1608 write (iout,*) "Finish INT_TO_CART"
1611 dc(j,i)=c(j,i+1)-c(j,i)
1612 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1616 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
1618 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1619 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1625 write(iout,*) "read angles from input"
1626 call read_angles(inp,*36)
1631 36 write (iout,'(a)') 'Error reading angle file.'
1633 call mpi_finalize( MPI_COMM_WORLD,IERR )
1635 stop 'Error reading angle file.'
1637 else if (extconf) then
1638 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0) &
1639 write (iout,'(a)') 'Extended chain initial geometry.'
1641 theta(i)=90d0*deg2rad
1642 if (molnum(i).eq.2) theta(i)=160d0*deg2rad
1645 phi(i)=180d0*deg2rad
1646 if (molnum(i).eq.2) phi(i)=30d0*deg2rad
1649 alph(i)=110d0*deg2rad
1650 if (molnum(i).eq.2) alph(i)=30d0*deg2rad
1653 omeg(i)=-120d0*deg2rad
1654 if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
1655 if (itype(i,1).le.0) omeg(i)=-omeg(i)
1659 if(me.eq.king.or..not.out1file) &
1660 write (iout,'(a)') 'Random-generated initial geometry.'
1664 if (me.eq.king .or. fg_rank.eq.0 .and. &
1665 ( modecalc.eq.12 .or. modecalc.eq.14) ) then
1669 call gen_rand_conf(itmp,*30)
1671 30 write (iout,*) 'Failed to generate random conformation',&
1673 write (*,*) 'Processor:',me,&
1674 ' Failed to generate random conformation',&
1684 write (iout,'(a,i3,a)') 'Processor:',me,&
1685 ' error in generating random conformation.'
1686 write (*,'(a,i3,a)') 'Processor:',me,&
1687 ' error in generating random conformation.'
1690 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1696 call gen_rand_conf(itmp,*335)
1698 335 write (iout,*) 'Failed to generate random conformation',&
1700 write (*,*) 'Failed to generate random conformation',&
1703 write (iout,'(a,i3,a)') 'Processor:',me,&
1704 ' error in generating random conformation.'
1705 write (*,'(a,i3,a)') 'Processor:',me,&
1706 ' error in generating random conformation.'
1711 elseif (modecalc.eq.4) then
1712 read (inp,'(a)') intinname
1713 open (intin,file=intinname,status='old',err=333)
1714 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0) &
1715 write (iout,'(a)') 'intinname',intinname
1716 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1718 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1720 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1722 stop 'Error opening angle file.'
1726 ! Generate distance constraints, if the PDB structure is to be regularized.
1727 if (nthread.gt.0) then
1728 call read_threadbase
1731 if (me.eq.king .or. .not. out1file) &
1733 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1734 write (iout,'(/a,i3,a)') &
1735 'The chain contains',ns,' disulfide-bridging cysteines.'
1736 write (iout,'(20i4)') (iss(i),i=1,ns)
1738 write(iout,*)"Running with dynamic disulfide-bond formation"
1740 write (iout,'(/a/)') 'Pre-formed links are:'
1746 if (me.eq.king.or..not.out1file) &
1747 write (iout,'(2a,i3,3a,i3,a,3f10.3)') &
1748 restyp(it1,1),'(',i1,') -- ',restyp(it2,1),'(',i2,')',dhpb(i),&
1754 if (ns.gt.0.and.dyn_ss) then
1758 forcon(i-nss)=forcon(i)
1765 dyn_ss_mask(iss(i))=.true.
1768 if (i2ndstr.gt.0) call secstrp2dihc
1769 if (indpdb.gt.0) then
1770 write(iout,*) "WCHODZE TU!!"
1771 call int_from_cart1(.true.)
1773 ! call geom_to_var(nvar,x)
1774 ! call etotal(energia(0))
1775 ! call enerprint(energia(0))
1776 ! call briefout(0,etot)
1778 !d write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1779 !d write (iout,'(a)') 'Variable list:'
1780 !d write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1782 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file)) &
1783 write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') &
1784 'Processor',myrank,': end reading molecular data.'
1787 end subroutine molread
1788 !-----------------------------------------------------------------------------
1789 subroutine read_constr_homology
1790 use control, only:init_int_table,homology_partition
1791 use MD_data, only:iset
1793 ! include 'DIMENSIONS'
1797 ! include 'COMMON.SETUP'
1798 ! include 'COMMON.CONTROL'
1799 ! include 'COMMON.HOMOLOGY'
1800 ! include 'COMMON.CHAIN'
1801 ! include 'COMMON.IOUNITS'
1802 ! include 'COMMON.MD'
1803 ! include 'COMMON.QRESTR'
1804 ! include 'COMMON.GEO'
1805 ! include 'COMMON.INTERACT'
1806 ! include 'COMMON.NAMES'
1807 ! include 'COMMON.VAR'
1810 ! double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
1812 ! common /przechowalnia/ odl_temp(maxres,maxres,max_template),
1813 ! & sigma_odl_temp(maxres,maxres,max_template)
1815 character*24 model_ki_dist, model_ki_angle
1816 character*500 controlcard
1817 integer :: ki,i,ii,j,k,l
1818 integer, dimension (:), allocatable :: ii_in_use
1819 integer :: i_tmp,idomain_tmp,&
1820 irec,ik,iistart,nres_temp
1823 logical :: liiflag,lfirst
1826 ! FP - Nov. 2014 Temporary specifications for new vars
1828 real(kind=8) :: rescore_tmp,x12,y12,z12,rescore2_tmp,&
1829 rescore3_tmp, dist_cut
1830 real(kind=8), dimension (:,:),allocatable :: rescore
1831 real(kind=8), dimension (:,:),allocatable :: rescore2
1832 real(kind=8), dimension (:,:),allocatable :: rescore3
1833 real(kind=8) :: distal
1834 character*24 tpl_k_rescore
1835 character*256 pdbfile
1837 ! -----------------------------------------------------------------
1838 ! Reading multiple PDB ref structures and calculation of retraints
1839 ! not using pre-computed ones stored in files model_ki_{dist,angle}
1841 ! -----------------------------------------------------------------
1844 ! Alternative: reading from input
1845 call card_concat(controlcard,.true.)
1846 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
1847 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
1848 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
1849 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
1850 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
1851 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
1852 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
1853 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
1854 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
1855 if(.not.read2sigma.and.start_from_model) then
1856 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)&
1857 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
1858 start_from_model=.false.
1859 iranconf=(indpdb.le.0)
1861 if(start_from_model .and. (me.eq.king .or. .not. out1file))&
1862 write(iout,*) 'START_FROM_MODELS is ON'
1863 ! if(start_from_model .and. rest) then
1864 ! if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1865 ! write(iout,*) 'START_FROM_MODELS is OFF'
1866 ! write(iout,*) 'remove restart keyword from input'
1869 if (start_from_model) nmodel_start=constr_homology
1870 if(.not.allocated(waga_homology)) allocate (waga_homology(homol_nset))
1871 if (homol_nset.gt.1)then
1872 call card_concat(controlcard,.true.)
1873 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
1874 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
1875 ! write(iout,*) "iset homology_weight "
1877 write(iout,*) i,waga_homology(i)
1880 iset=mod(kolor,homol_nset)+1
1883 waga_homology(1)=1.0
1886 !d write (iout,*) "nnt",nnt," nct",nct
1889 if (read_homol_frag) then
1890 call read_klapaucjusz
1896 ! write(iout,*) 'nnt=',nnt,'nct=',nct
1899 ! do k=1,constr_homology
1913 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
1914 if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
1916 do k=1,constr_homology
1918 read(inp,'(a)') pdbfile
1919 pdbfiles_chomo(k)=pdbfile
1920 if(me.eq.king .or. .not. out1file) &
1921 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',&
1922 pdbfile(:ilen(pdbfile))
1923 open(ipdbin,file=pdbfile,status='old',err=33)
1925 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
1926 pdbfile(:ilen(pdbfile))
1929 ! print *,'Begin reading pdb data'
1931 ! Files containing res sim or local scores (former containing sigmas)
1934 write(kic2,'(bz,i2.2)') k
1936 tpl_k_rescore="template"//kic2//".sco"
1937 write(iout,*) "tpl_k_rescore=",tpl_k_rescore
1940 write(iout,*) "read2sigma",read2sigma
1942 if (read2sigma) then
1943 call readpdb_template(k)
1947 write(iout,*) "after readpdb"
1948 if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
1951 if(.not.allocated(rescore)) allocate(rescore(constr_homology,nres))
1952 if(.not.allocated(rescore2)) allocate(rescore2(constr_homology,nres))
1953 if(.not.allocated(rescore3)) allocate(rescore3(constr_homology,nres))
1954 if(.not.allocated(ii_in_use)) allocate(ii_in_use(nres*(nres-1)))
1955 if(.not.allocated(idomain)) allocate(idomain(constr_homology,nres))
1956 if(.not.allocated(l_homo)) allocate(l_homo(constr_homology,1000*nres))
1957 if(.not.allocated(ires_homo)) allocate(ires_homo(nres*200))
1958 if(.not.allocated(jres_homo)) allocate(jres_homo(nres*200))
1959 if(.not.allocated(odl)) allocate(odl(constr_homology,nres*200))
1960 if(.not.allocated(sigma_odl)) allocate(sigma_odl(constr_homology,nres*200))
1961 if(.not.allocated(dih)) allocate(dih(constr_homology,nres))
1962 if(.not.allocated(sigma_dih)) allocate(sigma_dih(constr_homology,nres))
1963 if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1964 if(.not.allocated(sigma_theta)) allocate(sigma_theta(constr_homology,nres))
1965 ! if(.not.allocated(thetatpl)) allocate(thetatpl(constr_homology,nres))
1966 if(.not.allocated(sigma_d)) allocate(sigma_d(constr_homology,nres))
1967 if(.not.allocated(xxtpl)) allocate(xxtpl(constr_homology,nres))
1968 if(.not.allocated(yytpl)) allocate(yytpl(constr_homology,nres))
1969 if(.not.allocated(zztpl)) allocate(zztpl(constr_homology,nres))
1970 ! if(.not.allocated(distance)) allocate(distance(constr_homology))
1971 ! if(.not.allocated(distancek)) allocate(distancek(constr_homology))
1975 ! Distance restraints
1978 ! Copy the coordinates from reference coordinates (?)
1979 do i=1,2*nres_chomo(k)
1982 ! write (iout,*) "c(",j,i,") =",c(j,i)
1986 ! From read_dist_constr (commented out 25/11/2014 <-> res sim)
1988 ! write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
1989 open (ientin,file=tpl_k_rescore,status='old')
1990 if (nnt.gt.1) rescore(k,1)=0.0d0
1991 do irec=nnt,nct ! loop for reading res sim
1992 if (read2sigma) then
1993 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,&
1994 rescore3_tmp,idomain_tmp
1996 idomain(k,i_tmp)=idomain_tmp
1997 rescore(k,i_tmp)=rescore_tmp
1998 rescore2(k,i_tmp)=rescore2_tmp
1999 rescore3(k,i_tmp)=rescore3_tmp
2000 if (.not. out1file .or. me.eq.king)&
2001 write(iout,'(a7,i5,3f10.5,i5)') "rescore",&
2002 i_tmp,rescore2_tmp,rescore_tmp,&
2003 rescore3_tmp,idomain_tmp
2006 read (ientin,*,end=1401) rescore_tmp
2008 ! rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2009 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2010 ! write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2015 if (waga_dist.ne.0.0d0) then
2023 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2024 ! write (iout,*) k,i,j,distal,dist2_cut
2026 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2027 .and. distal.le.dist2_cut ) then
2033 ! write (iout,*) "k",k
2034 ! write (iout,*) "i",i," j",j," constr_homology",
2039 if (read2sigma) then
2042 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2044 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2045 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2046 sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2048 if (odl(k,ii).le.dist_cut) then
2049 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2052 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2053 dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2055 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2056 dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2060 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2063 ! l_homo(k,ii)=.false.
2069 ! write (iout,*) "Distance restraints set"
2072 ! Theta, dihedral and SC retraints
2074 if (waga_angle.gt.0.0d0) then
2075 ! open (ientin,file=tpl_k_sigma_dih,status='old')
2076 ! do irec=1,maxres-3 ! loop for reading sigma_dih
2077 ! read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2078 ! if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2079 ! sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2080 ! & sigma_dih(k,i+nnt-1)
2085 if (idomain(k,i).eq.0) then
2089 dih(k,i)=phiref(i) ! right?
2090 ! read (ientin,*) sigma_dih(k,i) ! original variant
2091 ! write (iout,*) "dih(",k,i,") =",dih(k,i)
2092 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2093 ! & "rescore(",k,i-1,") =",rescore(k,i-1),
2094 ! & "rescore(",k,i-2,") =",rescore(k,i-2),
2095 ! & "rescore(",k,i-3,") =",rescore(k,i-3)
2097 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2098 rescore(k,i-2)+rescore(k,i-3))/4.0
2099 ! if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2100 ! write (iout,*) "Raw sigmas for dihedral angle restraints"
2101 ! write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2102 ! sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2103 ! rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2104 ! Instead of res sim other local measure of b/b str reliability possible
2105 if (sigma_dih(k,i).ne.0) &
2106 sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2107 ! sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2111 ! write (iout,*) "Dihedral angle restraints set"
2114 if (waga_theta.gt.0.0d0) then
2115 ! open (ientin,file=tpl_k_sigma_theta,status='old')
2116 ! do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2117 ! read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2118 ! sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2119 ! & sigma_theta(k,i+nnt-1)
2124 do i = nnt+2,nct ! right? without parallel.
2125 ! do i = i=1,nres ! alternative for bounds acc to readpdb?
2126 ! do i=ithet_start,ithet_end ! with FG parallel.
2127 if (idomain(k,i).eq.0) then
2128 sigma_theta(k,i)=0.0
2131 thetatpl(k,i)=thetaref(i)
2132 ! write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2133 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2134 ! & "rescore(",k,i-1,") =",rescore(k,i-1),
2135 ! & "rescore(",k,i-2,") =",rescore(k,i-2)
2136 ! read (ientin,*) sigma_theta(k,i) ! 1st variant
2137 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ &
2139 ! if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2140 if (sigma_theta(k,i).ne.0) &
2141 sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2143 ! sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2144 ! rescore(k,i-2) ! right expression ?
2145 ! sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2148 ! write (iout,*) "Angle restraints set"
2151 if (waga_d.gt.0.0d0) then
2152 ! open (ientin,file=tpl_k_sigma_d,status='old')
2153 ! do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2154 ! read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2155 ! sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2156 ! & sigma_d(k,i+nnt-1)
2160 do i = nnt,nct ! right? without parallel.
2161 ! do i=2,nres-1 ! alternative for bounds acc to readpdb?
2162 ! do i=loc_start,loc_end ! with FG parallel.
2163 if (itype(i,1).eq.10) cycle
2164 if (idomain(k,i).eq.0 ) then
2171 ! write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2172 ! write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2173 ! write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2174 ! write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2175 sigma_d(k,i)=rescore3(k,i) ! right expression ?
2176 if (sigma_d(k,i).ne.0) &
2177 sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2179 ! sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2180 ! sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2181 ! read (ientin,*) sigma_d(k,i) ! 1st variant
2185 ! write (iout,*) "SC restraints set"
2188 ! remove distance restraints not used in any model from the list
2189 ! shift data in all arrays
2191 ! write (iout,*) "waga_dist",waga_dist," nnt",nnt," nct",nct
2192 if (waga_dist.ne.0.0d0) then
2199 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2200 ! & .and. distal.le.dist2_cut ) then
2201 ! write (iout,*) "i",i," j",j," ii",ii
2203 if (ii_in_use(ii).eq.0.and.liiflag.or. &
2204 ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2211 if(i10.eq.lim_odl) i10=i10+1
2213 ires_homo(iistart+ki)=ires_homo(ki+i01)
2214 jres_homo(iistart+ki)=jres_homo(ki+i01)
2215 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2216 do k=1,constr_homology
2217 odl(k,iistart+ki)=odl(k,ki+i01)
2218 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2219 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2222 iistart=iistart+i10-i01
2225 if (ii_in_use(ii).ne.0.and..not.liiflag) then
2233 ! write (iout,*) "Removing distances completed"
2235 endif ! .not. klapaucjusz
2237 if (constr_homology.gt.0) call homology_partition
2238 write (iout,*) "After homology_partition"
2240 if (constr_homology.gt.0) call init_int_table
2241 write (iout,*) "After init_int_table"
2243 ! endif ! .not. klapaucjusz
2245 ! if (constr_homology.gt.0) call homology_partition
2246 ! write (iout,*) "After homology_partition"
2248 ! if (constr_homology.gt.0) call init_int_table
2249 ! write (iout,*) "After init_int_table"
2251 ! write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2252 ! write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2256 if (.not.out_template_restr) return
2257 !d write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2258 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2259 write (iout,*) "Distance restraints from templates"
2261 write(iout,'(3i7,100(2f8.2,1x,l1,4x))') &
2262 ii,ires_homo(ii),jres_homo(ii),&
2263 (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),&
2264 ki=1,constr_homology)
2266 write (iout,*) "Dihedral angle restraints from templates"
2268 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2269 (rad2deg*dih(ki,i),&
2270 rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2272 write (iout,*) "Virtual-bond angle restraints from templates"
2274 write (iout,'(i7,a4,100(2f8.2,4x))') i,restyp(itype(i,1),1),&
2275 (rad2deg*thetatpl(ki,i),&
2276 rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2278 write (iout,*) "SC restraints from templates"
2280 write(iout,'(i7,100(4f8.2,4x))') i,&
2281 (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), &
2282 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2286 end subroutine read_constr_homology
2287 !-----------------------------------------------------------------------------
2288 subroutine read_klapaucjusz
2291 ! include 'DIMENSIONS'
2295 ! include 'COMMON.SETUP'
2296 ! include 'COMMON.CONTROL'
2297 ! include 'COMMON.HOMOLOGY'
2298 ! include 'COMMON.CHAIN'
2299 ! include 'COMMON.IOUNITS'
2300 ! include 'COMMON.MD'
2301 ! include 'COMMON.GEO'
2302 ! include 'COMMON.INTERACT'
2303 ! include 'COMMON.NAMES'
2304 character(len=256) fragfile
2305 integer, dimension(:), allocatable :: ninclust,nresclust,itype_temp,&
2307 integer, dimension(:,:), allocatable :: iresclust,inclust
2310 character(len=2) :: kic2
2311 character(len=24) :: model_ki_dist, model_ki_angle
2312 character(len=500) :: controlcard
2313 integer :: ki, i, j, jj,k, l, i_tmp,&
2315 ik,ll,lll,ii_old,ii,iii,ichain,kk,iistart,iishift,lim_xx,igr,&
2316 i01,i10,nnt_chain,nct_chain
2317 real(kind=8) :: distal
2318 logical :: lprn = .true.
2319 integer :: nres_temp
2322 logical :: liiflag,lfirst
2324 real(kind=8) rescore_tmp,x12,y12,z12,rescore2_tmp,dist_cut
2325 real(kind=8), dimension (:,:), allocatable :: rescore
2326 real(kind=8), dimension (:,:), allocatable :: rescore2
2327 character(len=24) :: tpl_k_rescore
2328 character(len=256) :: pdbfile
2331 ! For new homol impl
2333 ! include 'COMMON.VAR'
2335 ! write (iout,*) "READ_KLAPAUCJUSZ"
2336 ! print *,"READ_KLAPAUCJUSZ"
2338 call getenv("FRAGFILE",fragfile)
2339 write (iout,*) "Opening", fragfile
2341 open(ientin,file=fragfile,status="old",err=10)
2342 ! write (iout,*) " opened"
2351 itype_temp(:)=itype(:,1)
2355 ! write (iout,*) "Entering loop"
2358 DO IGR = 1,NCHAIN_GROUP
2360 ! write (iout,*) "igr",igr
2362 read(ientin,*) constr_homology,nclust
2363 if (start_from_model) then
2364 nmodel_start=constr_homology
2371 ichain=iequiv(1,igr)
2372 nnt_chain=chain_border(1,ichain)-chain_border1(1,ichain)+1
2373 nct_chain=chain_border(2,ichain)-chain_border1(1,ichain)+1
2374 ! write (iout,*) "nnt_chain",nnt_chain," nct_chain",nct_chain
2377 if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
2378 if(.not.allocated(nres_chomo)) allocate(nres_chomo(constr_homology))
2379 do k=1,constr_homology
2380 read(ientin,'(a)') pdbfile
2381 write (iout,'(a,5x,a)') 'KLAPAUCJUSZ: Opening PDB file', &
2382 pdbfile(:ilen(pdbfile))
2383 pdbfiles_chomo(k)=pdbfile
2384 open(ipdbin,file=pdbfile,status='old',err=33)
2386 33 write (iout,'(a,5x,a)') 'Error opening PDB file',&
2387 pdbfile(:ilen(pdbfile))
2392 call readpdb_template(k)
2402 read(ientin,*) ninclust(i),nresclust(i)
2403 read(ientin,*) (inclust(k,i),k=1,ninclust(i))
2404 read(ientin,*) (iresclust(k,i),k=1,nresclust(i))
2407 ! Loop over clusters
2410 do ll = 1,ninclust(l)
2413 ! write (iout,*) "l",l," ll",ll," k",k
2419 idomain(k,iresclust(i,l)+1) = 1
2421 idomain(k,iresclust(i,l)) = 1
2425 ! Distance restraints
2428 ! Copy the coordinates from reference coordinates (?)
2434 ! write (iout,*) "c(",j,i,") =",c(j,i)
2437 call int_from_cart(.true.,.false.)
2438 call sc_loc_geom(.false.)
2440 thetaref(i)=theta(i)
2444 if (waga_dist.ne.0.0d0) then
2447 do i = nnt_chain,nct_chain-2
2454 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2455 ! write (iout,*) k,i,j,distal,dist2_cut
2457 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 &
2458 .and. distal.le.dist2_cut ) then
2464 ! write (iout,*) "k",k
2465 ! write (iout,*) "i",i," j",j," constr_homology",
2467 ires_homo(ii)=i+chain_border1(1,igr)-1
2468 jres_homo(ii)=j+chain_border1(1,igr)-1
2470 if (read2sigma) then
2473 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2475 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2476 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = &
2477 sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2479 if (odl(k,ii).le.dist_cut) then
2480 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2483 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2484 dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2486 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* &
2487 dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2491 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2494 ! l_homo(k,ii)=.false.
2501 ! Theta, dihedral and SC retraints
2503 if (waga_angle.gt.0.0d0) then
2504 do i = nnt_chain+3,nct_chain
2505 iii=i+chain_border1(1,igr)-1
2506 if (idomain(k,i).eq.0) then
2507 ! sigma_dih(k,i)=0.0
2510 dih(k,iii)=phiref(i)
2512 (rescore(k,i)+rescore(k,i-1)+ &
2513 rescore(k,i-2)+rescore(k,i-3))/4.0
2514 ! write (iout,*) "k",k," l",l," i",i," rescore",rescore(k,i),
2515 ! & " sigma_dihed",sigma_dih(k,i)
2516 if (sigma_dih(k,iii).ne.0) &
2517 sigma_dih(k,iii)=1.0d0/(sigma_dih(k,iii)*sigma_dih(k,iii))
2522 if (waga_theta.gt.0.0d0) then
2523 do i = nnt_chain+2,nct_chain
2524 iii=i+chain_border1(1,igr)-1
2525 if (idomain(k,i).eq.0) then
2526 ! sigma_theta(k,i)=0.0
2529 thetatpl(k,iii)=thetaref(i)
2530 sigma_theta(k,iii)=(rescore(k,i)+rescore(k,i-1)+ &
2532 if (sigma_theta(k,iii).ne.0) &
2533 sigma_theta(k,iii)=1.0d0/ &
2534 (sigma_theta(k,iii)*sigma_theta(k,iii))
2538 if (waga_d.gt.0.0d0) then
2539 do i = nnt_chain,nct_chain
2540 iii=i+chain_border1(1,igr)-1
2541 if (itype(i,1).eq.10) cycle
2542 if (idomain(k,i).eq.0 ) then
2546 xxtpl(k,iii)=xxref(i)
2547 yytpl(k,iii)=yyref(i)
2548 zztpl(k,iii)=zzref(i)
2549 sigma_d(k,iii)=rescore(k,i)
2550 if (sigma_d(k,iii).ne.0) &
2551 sigma_d(k,iii)=1.0d0/(sigma_d(k,iii)*sigma_d(k,iii))
2552 ! if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1
2558 ! remove distance restraints not used in any model from the list
2559 ! shift data in all arrays
2561 ! write (iout,*) "ii_old",ii_old
2562 if (waga_dist.ne.0.0d0) then
2564 write (iout,*) "Distance restraints from templates"
2566 write(iout,'(4i5,100(2f8.2,1x,l1,4x))') &
2567 iii,ii_in_use(iii),ires_homo(iii),jres_homo(iii), &
2568 (odl(ki,iii),1.0d0/dsqrt(sigma_odl(ki,iii)),l_homo(ki,iii), &
2569 ki=1,constr_homology)
2575 do i=nnt_chain,nct_chain-2
2578 ! if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2579 ! & .and. distal.le.dist2_cut ) then
2580 ! write (iout,*) "i",i," j",j," ii",ii," i_in_use",ii_in_use(ii)
2582 if (ii_in_use(ii).eq.0.and.liiflag.or. &
2583 ii_in_use(ii).eq.1.and.liiflag.and.ii.eq.lim_odl) then
2590 if(i10.eq.lim_odl) i10=i10+1
2592 ires_homo(iistart+ki)=ires_homo(ki+i01)
2593 jres_homo(iistart+ki)=jres_homo(ki+i01)
2594 ii_in_use(iistart+ki)=ii_in_use(ki+i01)
2595 do k=1,constr_homology
2596 odl(k,iistart+ki)=odl(k,ki+i01)
2597 sigma_odl(k,iistart+ki)=sigma_odl(k,ki+i01)
2598 l_homo(k,iistart+ki)=l_homo(k,ki+i01)
2601 iistart=iistart+i10-i01
2604 if (ii_in_use(ii).ne.0.and..not.liiflag) then
2617 ichain=iequiv(i,igr)
2619 do j=nnt_chain,nct_chain
2620 jj=j+chain_border1(1,ichain)-chain_border1(1,iequiv(1,igr))
2621 do k=1,constr_homology
2623 sigma_dih(k,jj)=sigma_dih(k,j)
2624 thetatpl(k,jj)=thetatpl(k,j)
2625 sigma_theta(k,jj)=sigma_theta(k,j)
2626 xxtpl(k,jj)=xxtpl(k,j)
2627 yytpl(k,jj)=yytpl(k,j)
2628 zztpl(k,jj)=zztpl(k,j)
2629 sigma_d(k,jj)=sigma_d(k,j)
2633 jj=chain_border1(1,ichain)-chain_border1(1,iequiv(i-1,igr))
2634 ! write (iout,*) "igr",igr," i",i," ichain",ichain," jj",jj
2635 do j=ii_old+1,lim_odl
2636 ires_homo(j+lll)=ires_homo(j)+jj
2637 jres_homo(j+lll)=jres_homo(j)+jj
2638 do k=1,constr_homology
2639 odl(k,j+lll)=odl(k,j)
2640 sigma_odl(k,j+lll)=sigma_odl(k,j)
2641 l_homo(k,j+lll)=l_homo(k,j)
2652 if (waga_angle.gt.0.0d0) lim_dih=nct-nnt-2
2654 itype(:,1)=itype_temp(:)
2657 10 stop "Error in fragment file"
2658 end subroutine read_klapaucjusz
2660 !-----------------------------------------------------------
2661 subroutine seq2chains
2663 !c Split the total UNRES sequence, which has dummy residues separating
2664 !c the chains, into separate chains. The length of chain ichain is
2665 !c contained in chain_length(ichain), the first and last non-dummy
2666 !c residues are in chain_border(1,ichain) and chain_border(2,ichain),
2667 !c respectively. The lengths pertain to non-dummy residues only.
2670 ! include 'DIMENSIONS'
2671 use energy_data, only:molnum,nchain,chain_length,ireschain
2673 ! integer ireschain(nres)
2674 integer ii,ichain,i,j,mnum
2679 if (.not.allocated(chain_length)) allocate(chain_length(50))
2680 if (.not.allocated(chain_border)) allocate(chain_border(2,50))
2682 chain_length(ichain)=0
2684 do while (ii.lt.nres)
2685 write(iout,*) "in seq2chains",ii,new_chain
2687 if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
2688 if (.not.new_chain) then
2690 chain_border(2,ichain)=ii-1
2692 chain_border(1,ichain)=ii+1
2693 chain_length(ichain)=0
2697 chain_border(1,ichain)=ii
2700 chain_length(ichain)=chain_length(ichain)+1
2704 if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
2707 chain_length(ichain)=chain_length(ichain)+1
2709 if (chain_length(ichain).gt.0) then
2710 chain_border(2,ichain)=ii
2717 do j=chain_border(1,i),chain_border(2,i)
2723 !---------------------------------------------------------------------
2724 subroutine chain_symmetry(npermchain,tabpermchain)
2726 !c Determine chain symmetry. nperm is the number of permutations and
2727 !c tabperchain contains the allowed permutations of the chains.
2730 ! include "DIMENSIONS"
2731 ! include "COMMON.IOUNITS"
2734 npermchain,tabpermchain(50,5040),&
2735 tabperm(50,5040),mapchain(50),&
2736 iequiv(50,nres),iflag(nres)
2737 integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
2738 nperm,npermc,ind,mnum
2739 if (nchain.eq.1) then
2742 !c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
2746 !c Look for equivalent chains
2748 write(iout,*) "nchain",nchain
2750 write(iout,*) "chain",i," from",chain_border(1,i),&
2751 " to",chain_border(2,i)
2753 "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
2761 if (iflag(i).gt.0) cycle
2763 nchain_group=nchain_group+1
2765 iequiv(iieq,nchain_group)=i
2767 if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
2769 !c do while(k.lt.chain_length(i) .and.
2770 !c & itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
2771 do k=0,chain_length(i)-1
2774 if (itype(chain_border(1,i)+k,mnum).ne.&
2775 itype(chain_border(1,j)+k,mnum)) exit
2777 if (k.lt.chain_length(i)) cycle
2780 iequiv(iieq,nchain_group)=j
2782 nequiv(nchain_group)=iieq
2784 write(iout,*) "Number of equivalent chain groups:",nchain_group
2785 write(iout,*) "Equivalent chain groups"
2787 write(iout,*) "group",i," #members",nequiv(i)," chains",&
2788 (iequiv(j,i),j=1,nequiv(i))
2794 mapchain(ind)=iequiv(j,i)
2797 write (iout,*) "mapchain"
2799 write (iout,*) i,mapchain(i)
2803 call permut(nequiv(i),nperm,tabperm)
2809 tabpermchain(k,j)=iequiv(tabperm(k,j),i)
2814 npermchain=npermchain*nperm
2820 tabpermchain(l,ind)=tabpermchain(l,j)
2823 tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
2832 itemp(mapchain(j))=tabpermchain(j,i)
2835 tabpermchain(j,i)=itemp(j)
2838 write(iout,*) "Number of chain permutations",npermchain
2839 write(iout,*) "Permutations"
2841 write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
2845 !c---------------------------------------------------------------------
2846 integer function tperm(i,iperm,tabpermchain)
2848 ! include 'DIMENSIONS'
2850 integer tabpermchain(50,5040)
2854 tperm=tabpermchain(i,iperm)
2859 !-----------------------------------------------------------------------------