2 !-----------------------------------------------------------------------
6 !-----------------------------------------------------------------------------
7 ! Max. number of AA residues
8 integer,parameter :: maxres=6500!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 !-----------------------------------------------------------------------------
29 subroutine read_bridge
30 ! Read information about disulfide bridges.
31 use geometry_data, only: nres
33 use control_data, only:out1file
35 ! implicit real*8 (a-h,o-z)
36 ! include 'DIMENSIONS'
40 ! include 'COMMON.IOUNITS'
41 ! include 'COMMON.GEO'
42 ! include 'COMMON.VAR'
43 ! include 'COMMON.INTERACT'
44 ! include 'COMMON.LOCAL'
45 ! include 'COMMON.NAMES'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.FFIELD'
48 ! include 'COMMON.SBRIDGE'
49 ! include 'COMMON.HEADER'
50 ! include 'COMMON.CONTROL'
51 ! include 'COMMON.DBASE'
52 ! include 'COMMON.THREAD'
53 ! include 'COMMON.TIME1'
54 ! include 'COMMON.SETUP'
58 ! Read bridging residues.
64 read (inp,*) (iss(i),i=1,ns)
68 if(me.eq.king.or..not.out1file) &
69 write (iout,*) 'ns=',ns
71 write(iout,*) ' iss:',(iss(i),i=1,ns)
72 ! Check whether the specified bridging residues are cystines.
74 write(iout,*) i,iss(i)
75 if (itype(iss(i),1).ne.1) then
76 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') &
77 'Do you REALLY think that the residue ',&
78 restyp(itype(iss(i),1),1),i,&
79 ' can form a disulfide bridge?!!!'
80 write (*,'(2a,i3,a)') &
81 'Do you REALLY think that the residue ',&
82 restyp(itype(iss(i),1),1),i,&
83 ' can form a disulfide bridge?!!!'
85 call MPI_Finalize(MPI_COMM_WORLD,ierror)
91 if(.not.allocated(ihpb)) allocate(ihpb(ns))
92 if(.not.allocated(jhpb)) allocate(jhpb(ns))
94 ! Read preformed bridges.
98 if(.not.allocated(ihpb)) allocate(ihpb(nss))
99 if(.not.allocated(jhpb)) allocate(jhpb(nss))
100 read (inp,*) (ihpb(i),jhpb(i),i=1,nss)
102 if(.not.allocated(dhpb)) allocate(dhpb(nss))
103 if(.not.allocated(forcon)) allocate(forcon(nss))!(maxdim) !el maxdim=(maxres-1)*(maxres-2)/2
104 if(.not.allocated(dhpb1)) allocate(dhpb1(nss))
105 if(.not.allocated(ibecarb)) allocate(ibecarb(nss))
106 ! el Initialize the bridge array
111 !--------------------
113 write(iout,*)'nss=',nss !el,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
115 write(iout,*)'ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
117 ! Check if the residues involved in bridges are in the specified list of
121 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j) &
122 .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
123 write (iout,'(a,i3,a)') 'Disulfide pair',i,&
124 ' contains residues present in other pairs.'
125 write (*,'(a,i3,a)') 'Disulfide pair',i,&
126 ' contains residues present in other pairs.'
128 call MPI_Finalize(MPI_COMM_WORLD,ierror)
134 if (ihpb(i).eq.iss(j)) goto 10
136 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
139 if (jhpb(i).eq.iss(j)) goto 20
141 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
152 write(iout,*) "end read_bridge"
154 end subroutine read_bridge
155 !-----------------------------------------------------------------------------
156 subroutine read_x(kanal,*)
160 use geometry, only:int_from_cart1
161 ! implicit real*8 (a-h,o-z)
162 ! include 'DIMENSIONS'
163 ! include 'COMMON.GEO'
164 ! include 'COMMON.VAR'
165 ! include 'COMMON.CHAIN'
166 ! include 'COMMON.IOUNITS'
167 ! include 'COMMON.CONTROL'
168 ! include 'COMMON.LOCAL'
169 ! include 'COMMON.INTERACT'
170 ! Read coordinates from input
173 integer :: l,k,j,i,kanal
175 read(kanal,'(8f10.5)',end=10,err=10) &
176 ((c(l,k),l=1,3),k=1,nres),&
177 ((c(l,k+nres),l=1,3),k=nnt,nct)
180 c(j,2*nres)=c(j,nres)
182 call int_from_cart1(.false.)
185 dc(j,i)=c(j,i+1)-c(j,i)
186 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
190 if (itype(i,1).ne.10 .and. itype(i,1).ne.ntyp1) then
192 dc(j,i+nres)=c(j,i+nres)-c(j,i)
193 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
200 end subroutine read_x
201 !-----------------------------------------------------------------------------
202 subroutine read_threadbase
207 ! implicit real*8 (a-h,o-z)
208 ! include 'DIMENSIONS'
209 ! include 'COMMON.IOUNITS'
210 ! include 'COMMON.GEO'
211 ! include 'COMMON.VAR'
212 ! include 'COMMON.INTERACT'
213 ! include 'COMMON.LOCAL'
214 ! include 'COMMON.NAMES'
215 ! include 'COMMON.CHAIN'
216 ! include 'COMMON.FFIELD'
217 ! include 'COMMON.SBRIDGE'
218 ! include 'COMMON.HEADER'
219 ! include 'COMMON.CONTROL'
220 ! include 'COMMON.DBASE'
221 ! include 'COMMON.THREAD'
222 ! include 'COMMON.TIME1'
227 ! Read pattern database for threading.
229 allocate(cart_base(3,maxres_base,nseq)) !(3,maxres_base,maxseq)
230 allocate(nres_base(3,nseq)) !(3,maxseq)
231 allocate(str_nam(nseq)) !(maxseq)
233 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),&
234 nres_base(2,i),nres_base(3,i)
235 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,&
237 ! write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
238 ! & nres_base(2,i),nres_base(3,i)
239 ! write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
243 if (weidis.eq.0.0D0) weidis=0.1D0
252 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
253 write (iout,'(a,i5)') 'nexcl: ',nexcl
254 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
256 end subroutine read_threadbase
257 !-----------------------------------------------------------------------------
259 !el subroutine read_angles(kanal,iscor,energ,iprot,*)
260 subroutine read_angles(kanal,*)
264 ! implicit real*8 (a-h,o-z)
265 ! include 'DIMENSIONS'
266 ! include 'DIMENSIONS.ZSCOPT'
267 ! include 'COMMON.INTERACT'
268 ! include 'COMMON.SBRIDGE'
269 ! include 'COMMON.GEO'
270 ! include 'COMMON.VAR'
271 ! include 'COMMON.CHAIN'
272 ! include 'COMMON.IOUNITS'
273 character(len=80) :: lineh
274 integer :: iscor,iprot,ic
275 real(kind=8) :: energ
277 subroutine read_angles(kanal,*)
283 ! implicit real*8 (a-h,o-z)
284 ! include 'DIMENSIONS'
285 ! include 'COMMON.GEO'
286 ! include 'COMMON.VAR'
287 ! include 'COMMON.CHAIN'
288 ! include 'COMMON.IOUNITS'
289 ! include 'COMMON.CONTROL'
291 ! Read angles from input
296 read(kanal,'(a80)',end=10,err=10) lineh
297 read(lineh(:5),*,err=8) ic
298 read(lineh(6:),*,err=8) energ
301 print *,'error, assuming e=1d10',lineh
305 read(lineh(18:),*,end=10,err=10) nss
307 read (lineh(20:),*,end=10,err=10) &
308 (IHPB(I),JHPB(I),I=1,NSS),iscor
310 read (lineh(20:),*,end=10,err=10) (IHPB(I),JHPB(I),I=1,8)
311 read (kanal,*,end=10,err=10) (IHPB(I),JHPB(I),&
314 ! print *,"energy",energ," iscor",iscor
316 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
317 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
318 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
319 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
322 ! 9/7/01 avoid 180 deg valence angle
323 if (theta(i).gt.179.99d0) theta(i)=179.99d0
325 theta(i)=deg2rad*theta(i)
326 phi(i)=deg2rad*phi(i)
327 alph(i)=deg2rad*alph(i)
328 omeg(i)=deg2rad*omeg(i)
332 end subroutine read_angles
333 !-----------------------------------------------------------------------------
334 subroutine reada(rekord,lancuch,wartosc,default)
337 character*(*) :: rekord,lancuch
338 real(kind=8) :: wartosc,default
339 integer :: iread !,ilen
341 iread=index(rekord,lancuch)
346 iread=iread+ilen(lancuch)+1
347 read (rekord(iread:),*,err=10,end=10) wartosc
352 !-----------------------------------------------------------------------------
353 subroutine readi(rekord,lancuch,wartosc,default)
356 character*(*) :: rekord,lancuch
357 integer :: wartosc,default
358 integer :: iread !,ilen
360 iread=index(rekord,lancuch)
365 iread=iread+ilen(lancuch)+1
366 read (rekord(iread:),*,err=10,end=10) wartosc
371 !-----------------------------------------------------------------------------
372 subroutine multreadi(rekord,lancuch,tablica,dim,default)
376 integer :: tablica(dim),default
377 character*(*) :: rekord,lancuch
378 character(len=80) :: aux
379 integer :: iread !,ilen
384 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
385 if (iread.eq.0) return
386 iread=iread+ilen(lancuch)+1
387 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
389 end subroutine multreadi
390 !-----------------------------------------------------------------------------
391 subroutine multreada(rekord,lancuch,tablica,dim,default)
395 real(kind=8) :: tablica(dim),default
396 character*(*) :: rekord,lancuch
397 character(len=80) :: aux
398 integer :: iread !,ilen
403 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
404 if (iread.eq.0) return
405 iread=iread+ilen(lancuch)+1
406 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
408 end subroutine multreada
409 !-----------------------------------------------------------------------------
410 subroutine card_concat(card,to_upper)
412 ! dla UNRESA to_upper jest zawsze .true.
413 ! implicit real*8 (a-h,o-z)
414 ! include 'DIMENSIONS'
415 ! include 'COMMON.IOUNITS'
417 character(len=80) :: karta !,ucase
420 read (inp,'(a)') karta
421 if (to_upper) karta=ucase(karta)
423 do while (karta(80:80).eq.'&')
424 card=card(:ilen(card)+1)//karta(:79)
425 read (inp,'(a)') karta
426 if (to_upper) karta=ucase(karta)
428 card=card(:ilen(card)+1)//karta
430 end subroutine card_concat
431 !----------------------------------------------------------------------------
432 subroutine read_afminp
435 use control_data, only:out1file
437 character*320 afmcard
440 call card_concat(afmcard,.true.)
441 call readi(afmcard,"BEG",afmbeg,0)
442 call readi(afmcard,"END",afmend,0)
443 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
444 call reada(afmcard,"VEL",velAFMconst,0.0d0)
445 print *,'FORCE=' ,forceAFMconst
446 !------ NOW PROPERTIES FOR AFM
449 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
451 distafminit=dsqrt(distafminit)
452 print *,'initdist',distafminit
454 end subroutine read_afminp
455 !-----------------------------------------------------------------------------
456 subroutine read_dist_constr
459 use geometry, only: dist
463 ! implicit real*8 (a-h,o-z)
464 ! include 'DIMENSIONS'
468 ! include 'COMMON.SETUP'
469 ! include 'COMMON.CONTROL'
470 ! include 'COMMON.CHAIN'
471 ! include 'COMMON.IOUNITS'
472 ! include 'COMMON.SBRIDGE'
473 integer,dimension(2,100) :: ifrag_,ipair_
474 real(kind=8),dimension(100) :: wfrag_,wpair_
475 character(len=640) :: controlcard
478 integer :: i,k,j,ddjk,ii,jj,itemp
479 integer :: nfrag_,npair_,ndist_
480 real(kind=8) :: dist_cut
482 ! write (iout,*) "Calling read_dist_constr"
483 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
485 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
486 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
487 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
488 if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
489 if(.not.allocated(forcon)) allocate(forcon(maxdim))
490 if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
491 if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim))
492 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
493 call gen_dist_constr2
496 call card_concat(controlcard,.true.)
497 call readi(controlcard,"NFRAG",nfrag_,0)
498 call readi(controlcard,"NPAIR",npair_,0)
499 call readi(controlcard,"NDIST",ndist_,0)
500 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
501 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
502 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
503 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
504 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
505 ! write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
506 ! write (iout,*) "IFRAG"
508 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
510 ! write (iout,*) "IPAIR"
512 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
514 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
515 ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
516 ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
517 ! if(.not.allocated(forcon)) allocate(forcon(maxdim))
521 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
522 if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
523 ifrag_(2,i)=nstart_sup+nsup-1
524 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
526 if (wfrag_(i).gt.0.0d0) then
527 do j=ifrag_(1,i),ifrag_(2,i)-1
529 ! write (iout,*) "j",j," k",k
531 if (constr_dist.eq.1) then
536 forcon(nhpb)=wfrag_(i)
537 else if (constr_dist.eq.2) then
538 if (ddjk.le.dist_cut) then
543 forcon(nhpb)=wfrag_(i)
550 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
553 if (.not.out1file .or. me.eq.king) &
554 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
555 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
557 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
558 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
565 if (wpair_(i).gt.0.0d0) then
573 do j=ifrag_(1,ii),ifrag_(2,ii)
574 do k=ifrag_(1,jj),ifrag_(2,jj)
578 forcon(nhpb)=wpair_(i)
581 if (.not.out1file .or. me.eq.king) &
582 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
583 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
585 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
586 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
593 ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
594 ! if (forcon(nhpb+1).gt.0.0d0) then
596 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
597 if (constr_dist.eq.11) then
598 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
599 ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
600 !EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
601 fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0)
602 forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0)
605 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
606 ibecarb(i),forcon(nhpb+1)
608 if (forcon(nhpb+1).gt.0.0d0) then
610 if (ibecarb(i).gt.0) then
614 if (dhpb(nhpb).eq.0.0d0) &
615 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
619 if (.not.out1file .or. me.eq.king) &
620 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
621 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
623 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
624 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
630 end subroutine read_dist_constr
631 !-----------------------------------------------------------------------------
632 subroutine gen_dist_constr2
635 use geometry, only: dist
640 real(kind=8) :: distance
641 if (constr_dist.eq.11) then
642 do i=nstart_sup,nstart_sup+nsup-1
643 do j=i+2,nstart_sup+nsup-1
645 if (distance.le.15.0) then
647 ihpb(nhpb)=i+nstart_seq-nstart_sup
648 jhpb(nhpb)=j+nstart_seq-nstart_sup
649 forcon(nhpb)=sqrt(0.04*distance)
650 fordepth(nhpb)=sqrt(40.0/distance)
651 dhpb(nhpb)=distance-0.1d0
652 dhpb1(nhpb)=distance+0.1d0
655 if (.not.out1file .or. me.eq.king) &
656 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
657 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
659 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
660 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
666 do i=nstart_sup,nstart_sup+nsup-1
667 do j=i+2,nstart_sup+nsup-1
669 ihpb(nhpb)=i+nstart_seq-nstart_sup
670 jhpb(nhpb)=j+nstart_seq-nstart_sup
677 end subroutine gen_dist_constr2
679 !-----------------------------------------------------------------------------
691 !-----------------------------------------------------------------------------
692 subroutine copy_to_tmp(source)
694 ! include "DIMENSIONS"
695 ! include "COMMON.IOUNITS"
696 character*(*) :: source
697 character(len=256) :: tmpfile
701 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
702 inquire(file=tmpfile,exist=ex)
704 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
705 " to temporary directory..."
706 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
707 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
710 end subroutine copy_to_tmp
711 !-----------------------------------------------------------------------------
712 subroutine move_from_tmp(source)
714 ! include "DIMENSIONS"
715 ! include "COMMON.IOUNITS"
716 character*(*) :: source
719 write (*,*) "Moving ",source(:ilen(source)),&
720 " from temporary directory to working directory"
721 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
722 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
724 end subroutine move_from_tmp
725 !-----------------------------------------------------------------------------
727 !-----------------------------------------------------------------------------
728 ! $Date: 1994/10/12 17:24:21 $
731 logical function find_arg(ipos,line,errflag)
733 integer, parameter :: maxlen = 80
734 character(len=80) :: line
735 character(len=1) :: empty=' ',equal='='
738 ! This function returns .TRUE., if an argument follows keyword keywd; if so
739 ! IPOS will point to the first non-blank character of the argument. Returns
740 ! .FALSE., if no argument follows the keyword; in this case IPOS points
741 ! to the first non-blank character of the next keyword.
743 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
747 if (line(ipos:ipos).eq.equal) then
750 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
753 if (ipos.gt.maxlen) errflag=.true.
759 end function find_arg
760 !-----------------------------------------------------------------------------
761 logical function find_group(iunit,jout,key1)
763 character*(*) :: key1
764 character(len=80) :: karta !,ucase
765 integer :: iunit,jout
772 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
773 read (iunit,'(a)',end=10) karta
775 write (jout,'(2a)') '> ',karta(1:78)
778 10 find_group=.false.
780 end function find_group
782 !-----------------------------------------------------------------------------
783 logical function iblnk(charc)
784 character(len=1) :: charc
787 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
790 !-----------------------------------------------------------------------------
792 integer function ilen(string)
793 character*(*) :: string
797 1 if ( ilen .gt. 0 ) then
798 if ( iblnk( string(ilen:ilen) ) ) then
805 !-----------------------------------------------------------------------------
806 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
807 integer :: nkey,i,ikey,narg
808 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
809 ! character(len=16) :: ucase
812 if (ucase(keywd).eq.keywdset(i,ikey)) then
818 ! No match to the allowed set of keywords if this point is reached.
821 end function in_keywd_set
822 !-----------------------------------------------------------------------------
823 character function lcase(string)
824 integer :: i, k, idiff
825 character*(*) :: string
826 character(len=1) :: c
827 character(len=40) :: chtmp
833 if (string(k+1:) .ne. ' ') then
837 idiff = ichar('a') - ichar('A')
841 if (lge(c,'A') .and. lle(c,'Z')) then
842 lcase(i:i) = char(ichar(c) + idiff)
847 !-----------------------------------------------------------------------------
848 logical function lcom(ipos,karta)
849 character(len=80) :: karta
850 character :: koment(2) = (/'!','#'/)
855 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
859 !-----------------------------------------------------------------------------
860 logical function lower_case(ch)
862 lower_case=(ch.ge.'a' .and. ch.le.'z')
864 end function lower_case
865 !-----------------------------------------------------------------------------
866 subroutine mykey(line,keywd,ipos,blankline,errflag)
868 ! This subroutine seeks a non-empty substring keywd in the string LINE.
869 ! The substring begins with the first character different from blank and
870 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
871 ! at the character left to the first blank or "=". When the subroutine is
872 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
873 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
874 ! only separators or the maximum length of the data line (80) has been reached.
875 ! The logical variable ERRFLAG is set at .TRUE. if the string
876 ! consists only from a "=".
877 integer, parameter :: maxlen=80
878 character(len=1) :: empty=' ',equal='=',comma=','
879 character*(*) :: keywd
880 character(len=80) :: line
881 logical :: blankline,errflag !EL,lcom
882 integer :: ipos,istart,iend
884 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
887 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
888 ! At this point the rest of the input line turned out to contain only blanks
889 ! or to be commented out.
895 ! Checks whether the current char is a separator.
896 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
897 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
901 ! Error flag set to .true., if the length of the keyword was found less than 1.
902 if (iend.lt.istart) then
906 keywd=line(istart:iend)
909 !-----------------------------------------------------------------------------
910 subroutine numstr(inum,numm)
911 character(len=10) :: huj='0123456789'
912 character*(*) :: numm
913 integer :: inumm,inum,inum1,inum2
918 numm(3:3)=huj(inum2+1:inum2+1)
922 numm(2:2)=huj(inum2+1:inum2+1)
926 numm(1:1)=huj(inum2+1:inum2+1)
928 end subroutine numstr
930 !-----------------------------------------------------------------------------
931 function ucase(string)
932 integer :: i, k, idiff
933 character(*) :: string
934 character(len=len(string)) :: ucase
935 character(len=1) :: c
936 character(len=40) :: chtmp
942 if (string(k+1:) .ne. ' ') then
946 idiff = ichar('a') - ichar('A')
950 if (lge(c,'a') .and. lle(c,'z')) then
951 ucase(i:i) = char(ichar(c) - idiff)
956 !-----------------------------------------------------------------------------
958 !-----------------------------------------------------------------------------
959 subroutine pdbout(etot,tytul,iunit)
961 use geometry_data, only: c,nres
966 ! implicit real*8 (a-h,o-z)
967 ! include 'DIMENSIONS'
968 ! include 'COMMON.CHAIN'
969 ! include 'COMMON.INTERACT'
970 ! include 'COMMON.NAMES'
971 ! include 'COMMON.IOUNITS'
972 ! include 'COMMON.HEADER'
973 ! include 'COMMON.SBRIDGE'
974 ! include 'COMMON.DISTFIT'
975 ! include 'COMMON.MD'
976 !el character(len=50) :: tytul
977 character*(*) :: tytul
978 character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
980 integer,dimension(maxres) :: ica !(maxres)
982 integer,dimension(nres) :: ica !(maxres)
986 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
991 if(.not.allocated(molnum)) allocate(molnum(maxres))
993 if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
995 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
998 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
999 !model write (iunit,'(a5,i6)') 'MODEL',1
1001 if (nhfrag.gt.0) then
1003 iti=itype(hfrag(1,j),1)
1004 itj=itype(hfrag(2,j),1)
1006 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
1008 restyp(iti,1),hfrag(1,j)-1,&
1009 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1011 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
1013 restyp(iti,1),hfrag(1,j)-1,&
1014 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1019 if (nbfrag.gt.0) then
1023 iti=itype(bfrag(1,j),1)
1024 itj=itype(bfrag(2,j)-1,1)
1026 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1028 restyp(iti,1),bfrag(1,j)-1,&
1029 restyp(itj,1),bfrag(2,j)-2,0
1031 if (bfrag(3,j).gt.bfrag(4,j)) then
1033 itk=itype(bfrag(3,j),1)
1034 itl=itype(bfrag(4,j)+1,1)
1036 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)') &
1038 restyp(itl,1),bfrag(4,j),&
1039 restyp(itk,1),bfrag(3,j)-1,-1,&
1040 "N",restyp(itk,1),bfrag(3,j)-1,&
1041 "O",restyp(iti,1),bfrag(1,j)-1
1045 itk=itype(bfrag(3,j),1)
1046 itl=itype(bfrag(4,j)-1,1)
1049 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)') &
1051 restyp(itk,1),bfrag(3,j)-1,&
1052 restyp(itl,1),bfrag(4,j)-2,1,&
1053 "N",restyp(itk,1),bfrag(3,j)-1,&
1054 "O",restyp(iti,1),bfrag(1,j)-1
1066 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1067 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1068 'CYS',jdssb(i)-nnt+1
1070 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1071 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1072 'CYS',jhpb(i)-nnt+1-nres
1081 iti=itype(i,molnum(i))
1083 if (molnum(i+1).eq.0) then
1084 iti1=ntyp1_molec(molnum(i))
1086 iti1=itype(i+1,molnum(i+1))
1088 if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1090 if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1092 if (iti.eq.ntyp1_molec(molnum(i))) then
1095 write (iunit,'(a)') 'TER'
1100 if (molnum(i).eq.1) then
1102 write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1103 ires,(c(j,i),j=1,3),vtot(i)
1104 elseif(molnum(i).eq.2) then
1105 if (istype(i).eq.0) istype(i)=1
1106 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1107 chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1109 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1110 ires,(c(j,i),j=1,3),vtot(i)
1112 if ((iti.ne.10).and.(molnum(i).ne.5)) then
1114 if (molnum(i).eq.1) then
1115 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1116 ires,(c(j,nres+i),j=1,3),&
1118 else if (molnum(i).eq.2) then
1119 if (istype(i).eq.0) istype(i)=1
1120 write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1121 chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1127 write (iunit,'(a)') 'TER'
1129 if (itype(i,1).eq.ntyp1) cycle
1130 if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1131 write (iunit,30) ica(i),ica(i+1)
1132 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1133 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1134 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1135 write (iunit,30) ica(i),ica(i)+1
1138 if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1139 write (iunit,30) ica(nct),ica(nct)+1
1143 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1145 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1148 write (iunit,'(a6)') 'ENDMDL'
1149 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1151 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1152 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1153 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1155 30 FORMAT ('CONECT',8I5)
1156 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1159 end subroutine pdbout
1161 !-----------------------------------------------------------------------------
1162 subroutine MOL2out(etot,tytul)
1163 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1165 use geometry_data, only: c
1168 ! implicit real*8 (a-h,o-z)
1169 ! include 'DIMENSIONS'
1170 ! include 'COMMON.CHAIN'
1171 ! include 'COMMON.INTERACT'
1172 ! include 'COMMON.NAMES'
1173 ! include 'COMMON.IOUNITS'
1174 ! include 'COMMON.HEADER'
1175 ! include 'COMMON.SBRIDGE'
1176 character(len=32) :: tytul,fd
1177 character(len=3) :: zahl
1178 character(len=6) :: res_num,pom !,ucase
1182 real(kind=8) :: etot
1186 #elif (defined CRAY)
1191 write (imol2,'(a)') '#'
1192 write (imol2,'(a)') &
1193 '# Creating user name: unres'
1194 write (imol2,'(2a)') '# Creation time: ',&
1196 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1197 write (imol2,'(a)') tytul
1198 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1199 write (imol2,'(a)') 'SMALL'
1200 write (imol2,'(a)') 'USER_CHARGES'
1201 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1203 write (zahl,'(i3)') i
1204 pom=ucase(restyp(itype(i,1),1))
1205 res_num = pom(:3)//zahl(2:)
1206 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1208 write (imol2,'(a)') '\@<TRIPOS>BOND'
1210 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1213 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1215 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1217 write (zahl,'(i3)') i
1218 pom = ucase(restyp(itype(i,1),1))
1219 res_num = pom(:3)//zahl(2:)
1220 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1222 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1223 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1225 end subroutine MOL2out
1226 !-----------------------------------------------------------------------------
1230 use energy_data, only: itype
1232 ! implicit real*8 (a-h,o-z)
1233 ! include 'DIMENSIONS'
1234 ! include 'COMMON.IOUNITS'
1235 ! include 'COMMON.CHAIN'
1236 ! include 'COMMON.VAR'
1237 ! include 'COMMON.LOCAL'
1238 ! include 'COMMON.INTERACT'
1239 ! include 'COMMON.NAMES'
1240 ! include 'COMMON.GEO'
1241 ! include 'COMMON.TORSION'
1245 write (iout,'(/a)') 'Geometry of the virtual chain.'
1246 write (iout,'(7a)') ' Res ',' d',' Theta',&
1247 ' Phi',' Dsc',' Alpha',' Omega'
1250 ! print *,vbld(i),"vbld(i)",i
1251 write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
1252 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1256 end subroutine intout
1257 !-----------------------------------------------------------------------------
1259 subroutine briefout(it,ener,free)!,plik)
1261 subroutine briefout(it,ener)
1266 ! implicit real*8 (a-h,o-z)
1267 ! include 'DIMENSIONS'
1268 ! include 'COMMON.IOUNITS'
1269 ! include 'COMMON.CHAIN'
1270 ! include 'COMMON.VAR'
1271 ! include 'COMMON.LOCAL'
1272 ! include 'COMMON.INTERACT'
1273 ! include 'COMMON.NAMES'
1274 ! include 'COMMON.GEO'
1275 ! include 'COMMON.SBRIDGE'
1276 ! print '(a,i5)',intname,igeom
1279 real(kind=8) :: ener,free
1280 ! character(len=80) :: plik
1283 #if defined(AIX) || defined(PGI)
1284 open (igeom,file=intname,position='append')
1286 open (igeom,file=intname,access='append')
1295 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1297 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1299 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1301 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1303 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1305 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1306 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1307 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1308 ! if (nvar.gt.nphi+ntheta) then
1309 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1310 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1313 180 format (I5,F12.3,I2,9(1X,2I3))
1314 190 format (3X,11(1X,2I3))
1317 end subroutine briefout
1318 !-----------------------------------------------------------------------------
1320 subroutine fdate(fd)
1321 character(len=32) :: fd
1324 end subroutine fdate
1326 !-----------------------------------------------------------------------------
1328 real(kind=8) function gyrate(jcon)
1330 real(kind=8) function gyrate()
1333 use geometry_data, only: c
1336 ! implicit real*8 (a-h,o-z)
1337 ! include 'DIMENSIONS'
1338 ! include 'COMMON.INTERACT'
1339 ! include 'COMMON.CHAIN'
1341 real(kind=8),dimension(3) :: cen
1351 cen(j)=cen(j)+c(j,i)
1355 cen(j)=cen(j)/dble(nct-nnt+1)
1360 rg = rg + (c(j,i)-cen(j))**2
1364 gyrate = dsqrt(rg/dble(nct-nnt+1))
1366 gyrate = sqrt(rg/dble(nct-nnt+1))
1371 !-----------------------------------------------------------------------------
1373 subroutine reads(rekord,lancuch,wartosc,default)
1375 character*(*) :: rekord,lancuch,wartosc,default
1376 character(len=80) :: aux
1377 integer :: lenlan,lenrec,iread,ireade
1381 lenlan=ilen(lancuch)
1383 iread=index(rekord,lancuch(:lenlan)//"=")
1384 ! print *,"rekord",rekord," lancuch",lancuch
1385 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1386 if (iread.eq.0) then
1390 iread=iread+lenlan+1
1391 ! print *,"iread",iread
1392 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1393 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1395 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1397 ! print *,"iread",iread
1398 if (iread.gt.lenrec) then
1403 ! print *,"ireade",ireade
1404 do while (ireade.lt.lenrec .and. &
1405 .not.iblnk(rekord(ireade:ireade)))
1408 wartosc=rekord(iread:ireade)
1410 end subroutine reads
1412 !-----------------------------------------------------------------------------
1414 !-----------------------------------------------------------------------------
1415 subroutine permut(isym)
1417 use geometry_data, only: tabperm
1418 ! implicit real*8 (a-h,o-z)
1419 ! include 'DIMENSIONS'
1420 ! include 'COMMON.LOCAL'
1421 ! include 'COMMON.VAR'
1422 ! include 'COMMON.CHAIN'
1423 ! include 'COMMON.INTERACT'
1424 ! include 'COMMON.IOUNITS'
1425 ! include 'COMMON.GEO'
1426 ! include 'COMMON.NAMES'
1427 ! include 'COMMON.CONTROL'
1432 integer,dimension(isym) :: a
1433 ! parameter(n=symetr)
1446 10 print *,(a(i),i=1,n)
1450 ! write (iout,*) "tututu", kkk
1452 if(nextp(n,a)) go to 10
1454 end subroutine permut
1455 !-----------------------------------------------------------------------------
1456 logical function nextp(n,a)
1458 integer :: n,i,j,k,t
1460 integer,dimension(n) :: a
1462 10 if(a(i).lt.a(i+1)) go to 20
1479 if(a(j).lt.a(i)) go to 40
1487 !-----------------------------------------------------------------------------
1489 !-----------------------------------------------------------------------------
1490 integer function rescode(iseq,nam,itype,molecule)
1492 ! use io_base, only: ucase
1493 ! implicit real*8 (a-h,o-z)
1494 ! include 'DIMENSIONS'
1495 ! include 'COMMON.NAMES'
1496 ! include 'COMMON.IOUNITS'
1497 character(len=3) :: nam !,ucase
1498 integer :: iseq,itype,i
1500 print *,molecule,nam
1501 if (molecule.eq.1) then
1502 if (itype.eq.0) then
1504 do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1505 if (ucase(nam).eq.restyp(i,molecule)) then
1513 do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1514 if (nam(1:1).eq.onelet(i)) then
1521 else if (molecule.eq.2) then
1522 do i=1,ntyp1_molec(molecule)
1523 print *,nam(1:1),restyp(i,molecule)(1:1)
1524 if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
1529 else if (molecule.eq.3) then
1530 write(iout,*) "SUGAR not yet implemented"
1532 else if (molecule.eq.4) then
1533 write(iout,*) "Explicit LIPID not yet implemented"
1535 else if (molecule.eq.5) then
1536 do i=1,ntyp1_molec(molecule)
1537 print *,i,restyp(i,molecule)(1:2)
1538 if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then
1544 write(iout,*) "molecule not defined"
1546 write (iout,10) iseq,nam
1548 10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
1549 end function rescode
1550 integer function sugarcode(sugar,ires)
1553 if (sugar.eq.'D') then
1555 else if (sugar.eq.' ') then
1558 write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
1562 end function sugarcode
1563 !-----------------------------------------------------------------------------
1564 !-----------------------------------------------------------------------------