2 !-----------------------------------------------------------------------
6 !-----------------------------------------------------------------------------
7 ! Max. number of AA residues
8 integer,parameter :: maxres=6000!1200
9 ! Appr. max. number of interaction sites
10 integer,parameter :: maxres2=2*maxres
11 ! parameter (maxres6=6*maxres)
12 ! parameter (mmaxres2=(maxres2*(maxres2+1)/2))
13 !-----------------------------------------------------------------------------
14 ! Max. number of S-S bridges
15 ! integer,parameter :: maxss=20
16 !-----------------------------------------------------------------------------
17 ! Max. number of derivatives of virtual-bond and side-chain vectors in theta
19 ! integer,parameter :: maxdim=(maxres-1)*(maxres-2)/2
20 !-----------------------------------------------------------------------------
23 !-----------------------------------------------------------------------------
25 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
28 subroutine read_bridge
29 ! Read information about disulfide bridges.
30 use geometry_data, only: nres
32 use control_data, only:out1file
34 ! implicit real*8 (a-h,o-z)
35 ! include 'DIMENSIONS'
39 ! include 'COMMON.IOUNITS'
40 ! include 'COMMON.GEO'
41 ! include 'COMMON.VAR'
42 ! include 'COMMON.INTERACT'
43 ! include 'COMMON.LOCAL'
44 ! include 'COMMON.NAMES'
45 ! include 'COMMON.CHAIN'
46 ! include 'COMMON.FFIELD'
47 ! include 'COMMON.SBRIDGE'
48 ! include 'COMMON.HEADER'
49 ! include 'COMMON.CONTROL'
50 ! include 'COMMON.DBASE'
51 ! include 'COMMON.THREAD'
52 ! include 'COMMON.TIME1'
53 ! include 'COMMON.SETUP'
57 ! Read bridging residues.
63 read (inp,*) (iss(i),i=1,ns)
67 if(me.eq.king.or..not.out1file) &
68 write (iout,*) 'ns=',ns
70 write(iout,*) ' iss:',(iss(i),i=1,ns)
71 ! Check whether the specified bridging residues are cystines.
73 write(iout,*) i,iss(i)
74 if (itype(iss(i),1).ne.1) then
75 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
76 'Do you REALLY think that the residue ',&
77 restyp(itype(iss(i),1),1),i,&
78 ' can form a disulfide bridge?!!!'
79 write (*,'(2a,i3,a)') &
80 'Do you REALLY think that the residue ',&
81 restyp(itype(iss(i),1),1),i,&
82 ' can form a disulfide bridge?!!!'
84 call MPI_Finalize(MPI_COMM_WORLD,ierror)
90 if(.not.allocated(ihpb)) allocate(ihpb(ns))
91 if(.not.allocated(jhpb)) allocate(jhpb(ns))
93 ! Read preformed bridges.
97 if(.not.allocated(ihpb)) allocate(ihpb(nss))
98 if(.not.allocated(jhpb)) allocate(jhpb(nss))
99 read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
101 if(.not.allocated(dhpb)) allocate(dhpb(nss))
102 if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
103 if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
104 if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
105 ! el Initialize the bridge array
110 !--------------------
112 write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
114 write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
116 ! Check if the residues involved in bridges are in the specified list of
120 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
121 .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
122 write (iout,'(a,i3,a)') 'Disulfide pair',i,&
123 ' contains residues present in other pairs.'
124 write (*,'(a,i3,a)') 'Disulfide pair',i,&
125 ' contains residues present in other pairs.'
127 call MPI_Finalize(MPI_COMM_WORLD,ierror)
133 if (ihpb(i).eq.iss(j)) goto 10
135 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
138 if (jhpb(i).eq.iss(j)) goto 20
140 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
151 write(iout,*) "end read_bridge"
153 end subroutine read_bridge
154 !-----------------------------------------------------------------------------
155 subroutine read_x(kanal,*)
159 use geometry, only:int_from_cart1
160 ! implicit real*8 (a-h,o-z)
161 ! include 'DIMENSIONS'
162 ! include 'COMMON.GEO'
163 ! include 'COMMON.VAR'
164 ! include 'COMMON.CHAIN'
165 ! include 'COMMON.IOUNITS'
166 ! include 'COMMON.CONTROL'
167 ! include 'COMMON.LOCAL'
168 ! include 'COMMON.INTERACT'
169 ! Read coordinates from input
172 integer :: l,k,j,i,kanal
174 read(kanal,'(8f10.5)',end=10,err=10) &
175 ((c(l,k),l=1,3),k=1,nres),&
176 ((c(l,k+nres),l=1,3),k=nnt,nct)
179 c(j,2*nres)=c(j,nres)
181 call int_from_cart1(.false.)
184 dc(j,i)=c(j,i+1)-c(j,i)
185 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
189 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
191 dc(j,i+nres)=c(j,i+nres)-c(j,i)
192 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
199 end subroutine read_x
200 !-----------------------------------------------------------------------------
201 subroutine read_threadbase
206 ! implicit real*8 (a-h,o-z)
207 ! include 'DIMENSIONS'
208 ! include 'COMMON.IOUNITS'
209 ! include 'COMMON.GEO'
210 ! include 'COMMON.VAR'
211 ! include 'COMMON.INTERACT'
212 ! include 'COMMON.LOCAL'
213 ! include 'COMMON.NAMES'
214 ! include 'COMMON.CHAIN'
215 ! include 'COMMON.FFIELD'
216 ! include 'COMMON.SBRIDGE'
217 ! include 'COMMON.HEADER'
218 ! include 'COMMON.CONTROL'
219 ! include 'COMMON.DBASE'
220 ! include 'COMMON.THREAD'
221 ! include 'COMMON.TIME1'
226 ! Read pattern database for threading.
228 allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
229 allocate(nres_base(3,nseq)) !(3,maxseq)
230 allocate(str_nam(nseq)) !(maxseq)
232 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
233 nres_base(2,i),nres_base(3,i)
234 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
236 ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
237 ! & nres_base(2,i),nres_base(3,i)
238 ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
242 if (weidis.eq.0.0D0) weidis=0.1D0
251 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
252 write (iout,'(a,i5)') 'nexcl: ',nexcl
253 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
255 end subroutine read_threadbase
256 !-----------------------------------------------------------------------------
258 !el subroutine read_angles(kanal,iscor,energ,iprot,*)
259 subroutine read_angles(kanal,*)
263 ! implicit real*8 (a-h,o-z)
264 ! include 'DIMENSIONS'
265 ! include 'DIMENSIONS.ZSCOPT'
266 ! include 'COMMON.INTERACT'
267 ! include 'COMMON.SBRIDGE'
268 ! include 'COMMON.GEO'
269 ! include 'COMMON.VAR'
270 ! include 'COMMON.CHAIN'
271 ! include 'COMMON.IOUNITS'
272 character(len=80) :: lineh
273 integer :: iscor,iprot,ic
274 real(kind=8) :: energ
276 subroutine read_angles(kanal,*)
281 ! implicit real*8 (a-h,o-z)
282 ! include 'DIMENSIONS'
283 ! include 'COMMON.GEO'
284 ! include 'COMMON.VAR'
285 ! include 'COMMON.CHAIN'
286 ! include 'COMMON.IOUNITS'
287 ! include 'COMMON.CONTROL'
289 ! Read angles from input
294 read(kanal,'(a80)',end=10,err=10) lineh
295 read(lineh(:5),*,err=8) ic
296 read(lineh(6:),*,err=8) energ
299 print *,'error, assuming e=1d10',lineh
303 read(lineh(18:),*,end=10,err=10) nss
305 read (lineh(20:),*,end=10,err=10) &
306 (IHPB(I),JHPB(I),I=1,NSS),iscor
308 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
309 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
312 ! print *,"energy",energ," iscor",iscor
314 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
315 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
316 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
317 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
320 ! 9/7/01 avoid 180 deg valence angle
321 if (theta(i).gt.179.99d0) theta(i)=179.99d0
323 theta(i)=deg2rad*theta(i)
324 phi(i)=deg2rad*phi(i)
325 alph(i)=deg2rad*alph(i)
326 omeg(i)=deg2rad*omeg(i)
330 end subroutine read_angles
331 !-----------------------------------------------------------------------------
332 subroutine reada(rekord,lancuch,wartosc,default)
335 character*(*) :: rekord,lancuch
336 real(kind=8) :: wartosc,default
337 integer :: iread !,ilen
339 iread=index(rekord,lancuch)
344 iread=iread+ilen(lancuch)+1
345 read (rekord(iread:),*,err=10,end=10) wartosc
350 !-----------------------------------------------------------------------------
351 subroutine readi(rekord,lancuch,wartosc,default)
354 character*(*) :: rekord,lancuch
355 integer :: wartosc,default
356 integer :: iread !,ilen
358 iread=index(rekord,lancuch)
363 iread=iread+ilen(lancuch)+1
364 read (rekord(iread:),*,err=10,end=10) wartosc
369 !-----------------------------------------------------------------------------
370 subroutine multreadi(rekord,lancuch,tablica,dim,default)
374 integer :: tablica(dim),default
375 character*(*) :: rekord,lancuch
376 character(len=80) :: aux
377 integer :: iread !,ilen
382 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
383 if (iread.eq.0) return
384 iread=iread+ilen(lancuch)+1
385 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
387 end subroutine multreadi
388 !-----------------------------------------------------------------------------
389 subroutine multreada(rekord,lancuch,tablica,dim,default)
393 real(kind=8) :: tablica(dim),default
394 character*(*) :: rekord,lancuch
395 character(len=80) :: aux
396 integer :: iread !,ilen
401 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
402 if (iread.eq.0) return
403 iread=iread+ilen(lancuch)+1
404 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
406 end subroutine multreada
407 !-----------------------------------------------------------------------------
408 subroutine card_concat(card,to_upper)
410 ! dla UNRESA to_upper jest zawsze .true.
411 ! implicit real*8 (a-h,o-z)
412 ! include 'DIMENSIONS'
413 ! include 'COMMON.IOUNITS'
415 character(len=80) :: karta !,ucase
418 read (inp,'(a)') karta
419 if (to_upper) karta=ucase(karta)
421 do while (karta(80:80).eq.'&')
422 card=card(:ilen(card)+1)//karta(:79)
423 read (inp,'(a)') karta
424 if (to_upper) karta=ucase(karta)
426 card=card(:ilen(card)+1)//karta
428 end subroutine card_concat
429 !----------------------------------------------------------------------------
430 subroutine read_afminp
433 use control_data, only:out1file
435 character*320 afmcard
438 call card_concat(afmcard,.true.)
439 call readi(afmcard,"BEG",afmbeg,0)
440 call readi(afmcard,"END",afmend,0)
441 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
442 call reada(afmcard,"VEL",velAFMconst,0.0d0)
443 print *,'FORCE=' ,forceAFMconst
444 !------ NOW PROPERTIES FOR AFM
447 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
449 distafminit=dsqrt(distafminit)
450 print *,'initdist',distafminit
452 end subroutine read_afminp
453 !-----------------------------------------------------------------------------
454 subroutine read_dist_constr
457 use geometry, only: dist
461 ! implicit real*8 (a-h,o-z)
462 ! include 'DIMENSIONS'
466 ! include 'COMMON.SETUP'
467 ! include 'COMMON.CONTROL'
468 ! include 'COMMON.CHAIN'
469 ! include 'COMMON.IOUNITS'
470 ! include 'COMMON.SBRIDGE'
471 integer,dimension(2,100) :: ifrag_,ipair_
472 real(kind=8),dimension(100) :: wfrag_,wpair_
473 character(len=640) :: controlcard
476 integer :: i,k,j,ddjk,ii,jj,itemp
477 integer :: nfrag_,npair_,ndist_
478 real(kind=8) :: dist_cut
480 ! write (iout,*) "Calling read_dist_constr"
481 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
483 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
484 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
485 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
486 if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
487 if(.not.allocated(forcon)) allocate(forcon(maxdim))
488 if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
489 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
490 call gen_dist_constr2
493 call card_concat(controlcard,.true.)
494 call readi(controlcard,"NFRAG",nfrag_,0)
495 call readi(controlcard,"NPAIR",npair_,0)
496 call readi(controlcard,"NDIST",ndist_,0)
497 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
498 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
499 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
500 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
501 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
502 ! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
503 ! write (iout,*) "IFRAG"
505 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
507 ! write (iout,*) "IPAIR"
509 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
511 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
512 ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
513 ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
514 ! if(.not.allocated(forcon)) allocate(forcon(maxdim))
518 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
519 if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
520 ifrag_(2,i)=nstart_sup+nsup-1
521 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
523 if (wfrag_(i).gt.0.0d0) then
524 do j=ifrag_(1,i),ifrag_(2,i)-1
526 ! write (iout,*) "j",j," k",k
528 if (constr_dist.eq.1) then
533 forcon(nhpb)=wfrag_(i)
534 else if (constr_dist.eq.2) then
535 if (ddjk.le.dist_cut) then
540 forcon(nhpb)=wfrag_(i)
547 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
550 if (.not.out1file .or. me.eq.king) &
551 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
552 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
554 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
555 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
562 if (wpair_(i).gt.0.0d0) then
570 do j=ifrag_(1,ii),ifrag_(2,ii)
571 do k=ifrag_(1,jj),ifrag_(2,jj)
575 forcon(nhpb)=wpair_(i)
578 if (.not.out1file .or. me.eq.king) &
579 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
580 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
582 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
583 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
590 ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
591 ! if (forcon(nhpb+1).gt.0.0d0) then
593 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
594 if (constr_dist.eq.11) then
595 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
596 ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
597 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
600 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
601 ibecarb(i),forcon(nhpb+1)
603 if (forcon(nhpb+1).gt.0.0d0) then
605 if (ibecarb(i).gt.0) then
609 if (dhpb(nhpb).eq.0.0d0) &
610 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
614 if (.not.out1file .or. me.eq.king) &
615 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
616 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
618 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
619 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
625 end subroutine read_dist_constr
626 !-----------------------------------------------------------------------------
627 subroutine gen_dist_constr2
630 use geometry, only: dist
635 real(kind=8) :: distance
636 if (constr_dist.eq.11) then
637 do i=nstart_sup,nstart_sup+nsup-1
638 do j=i+2,nstart_sup+nsup-1
640 if (distance.le.15.0) then
642 ihpb(nhpb)=i+nstart_seq-nstart_sup
643 jhpb(nhpb)=j+nstart_seq-nstart_sup
644 forcon(nhpb)=sqrt(0.04*distance)
645 fordepth(nhpb)=sqrt(40.0/distance)
646 dhpb(nhpb)=distance-0.1d0
647 dhpb1(nhpb)=distance+0.1d0
650 if (.not.out1file .or. me.eq.king) &
651 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
652 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
654 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
655 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
661 do i=nstart_sup,nstart_sup+nsup-1
662 do j=i+2,nstart_sup+nsup-1
664 ihpb(nhpb)=i+nstart_seq-nstart_sup
665 jhpb(nhpb)=j+nstart_seq-nstart_sup
672 end subroutine gen_dist_constr2
674 !-----------------------------------------------------------------------------
686 !-----------------------------------------------------------------------------
687 subroutine copy_to_tmp(source)
689 ! include "DIMENSIONS"
690 ! include "COMMON.IOUNITS"
691 character*(*) :: source
692 character(len=256) :: tmpfile
696 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
697 inquire(file=tmpfile,exist=ex)
699 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
700 " to temporary directory..."
701 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
702 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
705 end subroutine copy_to_tmp
706 !-----------------------------------------------------------------------------
707 subroutine move_from_tmp(source)
709 ! include "DIMENSIONS"
710 ! include "COMMON.IOUNITS"
711 character*(*) :: source
714 write (*,*) "Moving ",source(:ilen(source)),&
715 " from temporary directory to working directory"
716 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
717 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
719 end subroutine move_from_tmp
720 !-----------------------------------------------------------------------------
722 !-----------------------------------------------------------------------------
723 ! $Date: 1994/10/12 17:24:21 $
726 logical function find_arg(ipos,line,errflag)
728 integer, parameter :: maxlen = 80
729 character(len=80) :: line
730 character(len=1) :: empty=' ',equal='='
733 ! This function returns .TRUE., if an argument follows keyword keywd; if so
734 ! IPOS will point to the first non-blank character of the argument. Returns
735 ! .FALSE., if no argument follows the keyword; in this case IPOS points
736 ! to the first non-blank character of the next keyword.
738 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
742 if (line(ipos:ipos).eq.equal) then
745 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
748 if (ipos.gt.maxlen) errflag=.true.
754 end function find_arg
755 !-----------------------------------------------------------------------------
756 logical function find_group(iunit,jout,key1)
758 character*(*) :: key1
759 character(len=80) :: karta !,ucase
760 integer :: iunit,jout
767 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
768 read (iunit,'(a)',end=10) karta
770 write (jout,'(2a)') '> ',karta(1:78)
773 10 find_group=.false.
775 end function find_group
776 !-----------------------------------------------------------------------------
777 logical function iblnk(charc)
778 character(len=1) :: charc
781 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
784 !-----------------------------------------------------------------------------
785 integer function ilen(string)
786 character*(*) :: string
790 1 if ( ilen .gt. 0 ) then
791 if ( iblnk( string(ilen:ilen) ) ) then
798 !-----------------------------------------------------------------------------
799 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
800 integer :: nkey,i,ikey,narg
801 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
802 ! character(len=16) :: ucase
805 if (ucase(keywd).eq.keywdset(i,ikey)) then
811 ! No match to the allowed set of keywords if this point is reached.
814 end function in_keywd_set
815 !-----------------------------------------------------------------------------
816 character function lcase(string)
817 integer :: i, k, idiff
818 character*(*) :: string
819 character(len=1) :: c
820 character(len=40) :: chtmp
826 if (string(k+1:) .ne. ' ') then
830 idiff = ichar('a') - ichar('A')
834 if (lge(c,'A') .and. lle(c,'Z')) then
835 lcase(i:i) = char(ichar(c) + idiff)
840 !-----------------------------------------------------------------------------
841 logical function lcom(ipos,karta)
842 character(len=80) :: karta
843 character :: koment(2) = (/'!','#'/)
848 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
852 !-----------------------------------------------------------------------------
853 logical function lower_case(ch)
855 lower_case=(ch.ge.'a' .and. ch.le.'z')
857 end function lower_case
858 !-----------------------------------------------------------------------------
859 subroutine mykey(line,keywd,ipos,blankline,errflag)
861 ! This subroutine seeks a non-empty substring keywd in the string LINE.
862 ! The substring begins with the first character different from blank and
863 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
864 ! at the character left to the first blank or "=". When the subroutine is
865 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
866 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
867 ! only separators or the maximum length of the data line (80) has been reached.
868 ! The logical variable ERRFLAG is set at .TRUE. if the string
869 ! consists only from a "=".
870 integer, parameter :: maxlen=80
871 character(len=1) :: empty=' ',equal='=',comma=','
872 character*(*) :: keywd
873 character(len=80) :: line
874 logical :: blankline,errflag !EL,lcom
875 integer :: ipos,istart,iend
877 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
880 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
881 ! At this point the rest of the input line turned out to contain only blanks
882 ! or to be commented out.
888 ! Checks whether the current char is a separator.
889 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
890 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
894 ! Error flag set to .true., if the length of the keyword was found less than 1.
895 if (iend.lt.istart) then
899 keywd=line(istart:iend)
902 !-----------------------------------------------------------------------------
903 subroutine numstr(inum,numm)
904 character(len=10) :: huj='0123456789'
905 character*(*) :: numm
906 integer :: inumm,inum,inum1,inum2
911 numm(3:3)=huj(inum2+1:inum2+1)
915 numm(2:2)=huj(inum2+1:inum2+1)
919 numm(1:1)=huj(inum2+1:inum2+1)
921 end subroutine numstr
922 !-----------------------------------------------------------------------------
923 function ucase(string)
924 integer :: i, k, idiff
925 character(*) :: string
926 character(len=len(string)) :: ucase
927 character(len=1) :: c
928 character(len=40) :: chtmp
934 if (string(k+1:) .ne. ' ') then
938 idiff = ichar('a') - ichar('A')
942 if (lge(c,'a') .and. lle(c,'z')) then
943 ucase(i:i) = char(ichar(c) - idiff)
948 !-----------------------------------------------------------------------------
950 !-----------------------------------------------------------------------------
951 subroutine pdbout(etot,tytul,iunit)
953 use geometry_data, only: c,nres
958 ! implicit real*8 (a-h,o-z)
959 ! include 'DIMENSIONS'
960 ! include 'COMMON.CHAIN'
961 ! include 'COMMON.INTERACT'
962 ! include 'COMMON.NAMES'
963 ! include 'COMMON.IOUNITS'
964 ! include 'COMMON.HEADER'
965 ! include 'COMMON.SBRIDGE'
966 ! include 'COMMON.DISTFIT'
967 ! include 'COMMON.MD'
968 !el character(len=50) :: tytul
969 character*(*) :: tytul
970 character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
971 integer,dimension(nres) :: ica !(maxres)
974 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
979 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
981 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
982 !model write (iunit,'(a5,i6)') 'MODEL',1
983 if (nhfrag.gt.0) then
985 iti=itype(hfrag(1,j),1)
986 itj=itype(hfrag(2,j),1)
988 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
990 restyp(iti,1),hfrag(1,j)-1,&
991 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
993 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
995 restyp(iti,1),hfrag(1,j)-1,&
996 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1001 if (nbfrag.gt.0) then
1005 iti=itype(bfrag(1,j),1)
1006 itj=itype(bfrag(2,j)-1,1)
1008 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1010 restyp(iti,1),bfrag(1,j)-1,&
1011 restyp(itj,1),bfrag(2,j)-2,0
1013 if (bfrag(3,j).gt.bfrag(4,j)) then
1015 itk=itype(bfrag(3,j),1)
1016 itl=itype(bfrag(4,j)+1,1)
1018 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
1020 restyp(itl,1),bfrag(4,j),&
1021 restyp(itk,1),bfrag(3,j)-1,-1,&
1022 "N",restyp(itk,1),bfrag(3,j)-1,&
1023 "O",restyp(iti,1),bfrag(1,j)-1
1027 itk=itype(bfrag(3,j),1)
1028 itl=itype(bfrag(4,j)-1,1)
1031 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3,2x,a1,2x,a3,i6,3x,a1,2x,a3,i6)') &
1033 restyp(itk,1),bfrag(3,j)-1,&
1034 restyp(itl,1),bfrag(4,j)-2,1,&
1035 "N",restyp(itk,1),bfrag(3,j)-1,&
1036 "O",restyp(iti,1),bfrag(1,j)-1
1048 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1049 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1050 'CYS',jdssb(i)-nnt+1
1052 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1053 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1054 'CYS',jhpb(i)-nnt+1-nres
1063 iti=itype(i,molnum(i))
1065 if (molnum(i+1).eq.0) then
1066 iti1=ntyp1_molec(molnum(i))
1068 iti1=itype(i+1,molnum(i+1))
1070 if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1072 if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1074 if (iti.eq.ntyp1_molec(molnum(i))) then
1077 write (iunit,'(a)') 'TER'
1082 if (molnum(i).eq.1) then
1084 write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1085 ires,(c(j,i),j=1,3),vtot(i)
1086 elseif(molnum(i).eq.2) then
1087 if (istype(i).eq.0) istype(i)=1
1088 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1089 chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1091 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1092 ires,(c(j,i),j=1,3),vtot(i)
1094 if ((iti.ne.10).and.(molnum(i).ne.5)) then
1096 if (molnum(i).eq.1) then
1097 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1098 ires,(c(j,nres+i),j=1,3),&
1100 else if (molnum(i).eq.2) then
1101 if (istype(i).eq.0) istype(i)=1
1102 write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1103 chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1109 write (iunit,'(a)') 'TER'
1111 if (itype(i,1).eq.ntyp1) cycle
1112 if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1113 write (iunit,30) ica(i),ica(i+1)
1114 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1115 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1116 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1117 write (iunit,30) ica(i),ica(i)+1
1120 if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1121 write (iunit,30) ica(nct),ica(nct)+1
1125 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1127 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1130 write (iunit,'(a6)') 'ENDMDL'
1131 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1133 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1134 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1135 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1137 30 FORMAT ('CONECT',8I5)
1138 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1141 end subroutine pdbout
1142 !-----------------------------------------------------------------------------
1143 subroutine MOL2out(etot,tytul)
1144 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1146 use geometry_data, only: c
1149 ! implicit real*8 (a-h,o-z)
1150 ! include 'DIMENSIONS'
1151 ! include 'COMMON.CHAIN'
1152 ! include 'COMMON.INTERACT'
1153 ! include 'COMMON.NAMES'
1154 ! include 'COMMON.IOUNITS'
1155 ! include 'COMMON.HEADER'
1156 ! include 'COMMON.SBRIDGE'
1157 character(len=32) :: tytul,fd
1158 character(len=3) :: zahl
1159 character(len=6) :: res_num,pom !,ucase
1163 real(kind=8) :: etot
1167 #elif (defined CRAY)
1172 write (imol2,'(a)') '#'
1173 write (imol2,'(a)') &
1174 '# Creating user name: unres'
1175 write (imol2,'(2a)') '# Creation time: ',&
1177 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1178 write (imol2,'(a)') tytul
1179 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1180 write (imol2,'(a)') 'SMALL'
1181 write (imol2,'(a)') 'USER_CHARGES'
1182 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1184 write (zahl,'(i3)') i
1185 pom=ucase(restyp(itype(i,1),1))
1186 res_num = pom(:3)//zahl(2:)
1187 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1189 write (imol2,'(a)') '\@<TRIPOS>BOND'
1191 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1194 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1196 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1198 write (zahl,'(i3)') i
1199 pom = ucase(restyp(itype(i,1),1))
1200 res_num = pom(:3)//zahl(2:)
1201 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1203 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1204 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1206 end subroutine MOL2out
1207 !-----------------------------------------------------------------------------
1211 use energy_data, only: itype
1213 ! implicit real*8 (a-h,o-z)
1214 ! include 'DIMENSIONS'
1215 ! include 'COMMON.IOUNITS'
1216 ! include 'COMMON.CHAIN'
1217 ! include 'COMMON.VAR'
1218 ! include 'COMMON.LOCAL'
1219 ! include 'COMMON.INTERACT'
1220 ! include 'COMMON.NAMES'
1221 ! include 'COMMON.GEO'
1222 ! include 'COMMON.TORSION'
1226 write (iout,'(/a)') 'Geometry of the virtual chain.'
1227 write (iout,'(7a)') ' Res ',' d',' Theta',&
1228 ' Phi',' Dsc',' Alpha',' Omega'
1231 ! print *,vbld(i),"vbld(i)",i
1232 write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1233 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1237 end subroutine intout
1238 !-----------------------------------------------------------------------------
1240 subroutine briefout(it,ener,free)!,plik)
1242 subroutine briefout(it,ener)
1247 ! implicit real*8 (a-h,o-z)
1248 ! include 'DIMENSIONS'
1249 ! include 'COMMON.IOUNITS'
1250 ! include 'COMMON.CHAIN'
1251 ! include 'COMMON.VAR'
1252 ! include 'COMMON.LOCAL'
1253 ! include 'COMMON.INTERACT'
1254 ! include 'COMMON.NAMES'
1255 ! include 'COMMON.GEO'
1256 ! include 'COMMON.SBRIDGE'
1257 ! print '(a,i5)',intname,igeom
1260 real(kind=8) :: ener,free
1261 ! character(len=80) :: plik
1264 #if defined(AIX) || defined(PGI)
1265 open (igeom,file=intname,position='append')
1267 open (igeom,file=intname,access='append')
1276 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1278 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1280 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1282 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1284 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1286 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1287 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1288 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1289 ! if (nvar.gt.nphi+ntheta) then
1290 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1291 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1294 180 format (I5,F12.3,I2,9(1X,2I3))
1295 190 format (3X,11(1X,2I3))
1298 end subroutine briefout
1299 !-----------------------------------------------------------------------------
1301 subroutine fdate(fd)
1302 character(len=32) :: fd
1305 end subroutine fdate
1307 !-----------------------------------------------------------------------------
1309 real(kind=8) function gyrate(jcon)
1311 real(kind=8) function gyrate()
1314 use geometry_data, only: c
1317 ! implicit real*8 (a-h,o-z)
1318 ! include 'DIMENSIONS'
1319 ! include 'COMMON.INTERACT'
1320 ! include 'COMMON.CHAIN'
1322 real(kind=8),dimension(3) :: cen
1332 cen(j)=cen(j)+c(j,i)
1336 cen(j)=cen(j)/dble(nct-nnt+1)
1341 rg = rg + (c(j,i)-cen(j))**2
1345 gyrate = dsqrt(rg/dble(nct-nnt+1))
1347 gyrate = sqrt(rg/dble(nct-nnt+1))
1352 !-----------------------------------------------------------------------------
1354 subroutine reads(rekord,lancuch,wartosc,default)
1356 character*(*) :: rekord,lancuch,wartosc,default
1357 character(len=80) :: aux
1358 integer :: lenlan,lenrec,iread,ireade
1362 lenlan=ilen(lancuch)
1364 iread=index(rekord,lancuch(:lenlan)//"=")
1365 ! print *,"rekord",rekord," lancuch",lancuch
1366 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1367 if (iread.eq.0) then
1371 iread=iread+lenlan+1
1372 ! print *,"iread",iread
1373 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1374 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1376 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1378 ! print *,"iread",iread
1379 if (iread.gt.lenrec) then
1384 ! print *,"ireade",ireade
1385 do while (ireade.lt.lenrec .and. &
1386 .not.iblnk(rekord(ireade:ireade)))
1389 wartosc=rekord(iread:ireade)
1391 end subroutine reads
1393 !-----------------------------------------------------------------------------
1395 !-----------------------------------------------------------------------------
1396 subroutine permut(isym)
1398 use geometry_data, only: tabperm
1399 ! implicit real*8 (a-h,o-z)
1400 ! include 'DIMENSIONS'
1401 ! include 'COMMON.LOCAL'
1402 ! include 'COMMON.VAR'
1403 ! include 'COMMON.CHAIN'
1404 ! include 'COMMON.INTERACT'
1405 ! include 'COMMON.IOUNITS'
1406 ! include 'COMMON.GEO'
1407 ! include 'COMMON.NAMES'
1408 ! include 'COMMON.CONTROL'
1413 integer,dimension(isym) :: a
1414 ! parameter(n=symetr)
1427 10 print *,(a(i),i=1,n)
1431 ! write (iout,*) "tututu", kkk
1433 if(nextp(n,a)) go to 10
1435 end subroutine permut
1436 !-----------------------------------------------------------------------------
1437 logical function nextp(n,a)
1439 integer :: n,i,j,k,t
1441 integer,dimension(n) :: a
1443 10 if(a(i).lt.a(i+1)) go to 20
1460 if(a(j).lt.a(i)) go to 40
1467 !-----------------------------------------------------------------------------
1468 !-----------------------------------------------------------------------------