2 !-----------------------------------------------------------------------
6 !-----------------------------------------------------------------------------
7 ! Max. number of AA residues
8 integer,parameter :: maxres=101000!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
438 real, dimension(3) ::cbeg,cend
441 call card_concat(afmcard,.true.)
442 call readi(afmcard,"BEG",afmbeg,0)
443 call readi(afmcard,"END",afmend,0)
444 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
445 call reada(afmcard,"VEL",velAFMconst,0.0d0)
446 print *,'FORCE=' ,forceAFMconst
449 if (afmbeg.eq.-1) then
450 read(inp,*) nbegafmmat,(afmbegcentr(i),i=1,nbegafmmat)
453 cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat
461 if (afmend.eq.-1) then
462 read(inp,*) nendafmmat,(afmendcentr(i),i=1,nendafmmat)
465 cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat
471 !------ NOW PROPERTIES FOR AFM
474 distafminit=(cend(i)-cbeg(i))**2+distafminit
476 distafminit=dsqrt(distafminit)
477 print *,'initdist',distafminit
479 end subroutine read_afminp
480 !C------------------------------------------------------------
481 subroutine read_afmnano
484 use control_data, only:out1file
488 read(inp,*,err=11) inanomove,(inanotab(i),i=1,inanomove),forcenanoconst
491 cm=cm+c(3,inanotab(i))
494 distnanoinit=cm-tubecenter(3)
498 ", number of center, their index and forceconstance"
501 end subroutine read_afmnano
503 !-----------------------------------------------------------------------------
504 subroutine read_dist_constr
507 use geometry, only: dist
511 ! implicit real*8 (a-h,o-z)
512 ! include 'DIMENSIONS'
516 ! include 'COMMON.SETUP'
517 ! include 'COMMON.CONTROL'
518 ! include 'COMMON.CHAIN'
519 ! include 'COMMON.IOUNITS'
520 ! include 'COMMON.SBRIDGE'
521 integer,dimension(2,100) :: ifrag_,ipair_
522 real(kind=8),dimension(100) :: wfrag_,wpair_
523 character(len=640) :: controlcard
526 integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1
527 integer :: nfrag_,npair_,ndist_
528 real(kind=8) :: dist_cut,ddjk
530 ! write (iout,*) "Calling read_dist_constr"
531 ! write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
533 if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
534 if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
535 if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
536 if(.not.allocated(dhpb1)) allocate(dhpb1(maxdim))
537 if(.not.allocated(forcon)) allocate(forcon(maxdim))
538 if(.not.allocated(fordepth)) allocate(fordepth(maxdim))
539 if(.not.allocated(ibecarb)) allocate(ibecarb(maxdim))
540 if ((genconstr.gt.0).and.(constr_dist.eq.11)) then
541 call gen_dist_constr2
544 call card_concat(controlcard,.true.)
545 call readi(controlcard,"NFRAG",nfrag_,0)
546 call readi(controlcard,"NPAIR",npair_,0)
547 call readi(controlcard,"NDIST",ndist_,0)
548 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
549 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
550 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
551 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
552 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
553 write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
554 ! write (iout,*) "IFRAG"
556 write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
558 ! write (iout,*) "IPAIR"
560 ! write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
562 ! if(.not.allocated(ihpb)) allocate(ihpb(maxdim))
563 ! if(.not.allocated(jhpb)) allocate(jhpb(maxdim))
564 ! if(.not.allocated(dhpb)) allocate(dhpb(maxdim))
565 ! if(.not.allocated(forcon)) allocate(forcon(maxdim))
569 ! if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
570 ! if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
571 ! ifrag_(2,i)=nstart_sup+nsup-1
572 ! write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
574 if (wfrag_(i).gt.0.0d0) then
575 do j=ifrag_(1,i),ifrag_(2,i)-1
577 write (iout,*) "j",j," k",k,nres
580 if (j.gt.nres) j1=j-nres
581 if (k.gt.nres) k1=k-nres
585 if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.&
586 (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle
587 write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj)
589 if (constr_dist.eq.1) then
594 forcon(nhpb)=wfrag_(i)
595 else if (constr_dist.eq.2) then
596 if (ddjk.le.dist_cut) then
601 forcon(nhpb)=wfrag_(i)
608 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
611 if (.not.out1file .or. me.eq.king) &
612 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
613 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
615 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",&
616 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
623 if (wpair_(i).gt.0.0d0) then
631 do j=ifrag_(1,ii),ifrag_(2,ii)
632 do k=ifrag_(1,jj),ifrag_(2,jj)
636 forcon(nhpb)=wpair_(i)
639 if (.not.out1file .or. me.eq.king) &
640 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
641 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
643 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
644 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
651 ! read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
652 ! if (forcon(nhpb+1).gt.0.0d0) then
654 ! dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
655 if (constr_dist.eq.11) then
656 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
657 ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
658 !EL fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
659 fordepth(nhpb+1)=fordepth(nhpb+1)**(0.25d0)
660 forcon(nhpb+1)=forcon(nhpb+1)**(0.25d0)
663 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), &
664 ibecarb(i),forcon(nhpb+1)
666 if (forcon(nhpb+1).gt.0.0d0) then
668 if (ibecarb(i).gt.0) then
672 if (dhpb(nhpb).eq.0.0d0) &
673 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
677 if (.not.out1file .or. me.eq.king) &
678 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
679 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
681 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",&
682 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
688 end subroutine read_dist_constr
689 !-----------------------------------------------------------------------------
690 subroutine gen_dist_constr2
693 use geometry, only: dist
698 real(kind=8) :: distance
699 if (constr_dist.eq.11) then
700 do i=nstart_sup,nstart_sup+nsup-1
701 do j=i+2,nstart_sup+nsup-1
703 if (distance.le.15.0) then
705 ihpb(nhpb)=i+nstart_seq-nstart_sup
706 jhpb(nhpb)=j+nstart_seq-nstart_sup
707 forcon(nhpb)=sqrt(0.04*distance)
708 fordepth(nhpb)=sqrt(40.0/distance)
709 dhpb(nhpb)=distance-0.1d0
710 dhpb1(nhpb)=distance+0.1d0
713 if (.not.out1file .or. me.eq.king) &
714 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
715 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
717 write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", &
718 nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
724 do i=nstart_sup,nstart_sup+nsup-1
725 do j=i+2,nstart_sup+nsup-1
727 ihpb(nhpb)=i+nstart_seq-nstart_sup
728 jhpb(nhpb)=j+nstart_seq-nstart_sup
735 end subroutine gen_dist_constr2
737 !-----------------------------------------------------------------------------
749 !-----------------------------------------------------------------------------
750 subroutine copy_to_tmp(source)
752 ! include "DIMENSIONS"
753 ! include "COMMON.IOUNITS"
754 character*(*) :: source
755 character(len=256) :: tmpfile
759 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
760 inquire(file=tmpfile,exist=ex)
762 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),&
763 " to temporary directory..."
764 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
765 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
768 end subroutine copy_to_tmp
769 !-----------------------------------------------------------------------------
770 subroutine move_from_tmp(source)
772 ! include "DIMENSIONS"
773 ! include "COMMON.IOUNITS"
774 character*(*) :: source
777 write (*,*) "Moving ",source(:ilen(source)),&
778 " from temporary directory to working directory"
779 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
780 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
782 end subroutine move_from_tmp
783 !-----------------------------------------------------------------------------
785 !-----------------------------------------------------------------------------
786 ! $Date: 1994/10/12 17:24:21 $
789 logical function find_arg(ipos,line,errflag)
791 integer, parameter :: maxlen = 80
792 character(len=80) :: line
793 character(len=1) :: empty=' ',equal='='
796 ! This function returns .TRUE., if an argument follows keyword keywd; if so
797 ! IPOS will point to the first non-blank character of the argument. Returns
798 ! .FALSE., if no argument follows the keyword; in this case IPOS points
799 ! to the first non-blank character of the next keyword.
801 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
805 if (line(ipos:ipos).eq.equal) then
808 do while (line(ipos:ipos) .eq. empty .and. ipos.le.maxlen)
811 if (ipos.gt.maxlen) errflag=.true.
817 end function find_arg
818 !-----------------------------------------------------------------------------
819 logical function find_group(iunit,jout,key1)
821 character*(*) :: key1
822 character(len=80) :: karta !,ucase
823 integer :: iunit,jout
830 do while (index(ucase(karta),key1(1:ll)).eq.0.or.lcom(1,karta))
831 read (iunit,'(a)',end=10) karta
833 write (jout,'(2a)') '> ',karta(1:78)
836 10 find_group=.false.
838 end function find_group
840 !-----------------------------------------------------------------------------
841 logical function iblnk(charc)
842 character(len=1) :: charc
845 iblnk = (n.eq.9) .or. (n.eq.10) .or. (charc.eq.' ')
848 !-----------------------------------------------------------------------------
850 integer function ilen(string)
851 character*(*) :: string
855 1 if ( ilen .gt. 0 ) then
856 if ( iblnk( string(ilen:ilen) ) ) then
863 !-----------------------------------------------------------------------------
864 integer function in_keywd_set(nkey,ikey,narg,keywd,keywdset)
865 integer :: nkey,i,ikey,narg
866 character(len=16) :: keywd,keywdset(1:nkey,0:nkey)
867 ! character(len=16) :: ucase
870 if (ucase(keywd).eq.keywdset(i,ikey)) then
876 ! No match to the allowed set of keywords if this point is reached.
879 end function in_keywd_set
880 !-----------------------------------------------------------------------------
881 character function lcase(string)
882 integer :: i, k, idiff
883 character*(*) :: string
884 character(len=1) :: c
885 character(len=40) :: chtmp
891 if (string(k+1:) .ne. ' ') then
895 idiff = ichar('a') - ichar('A')
899 if (lge(c,'A') .and. lle(c,'Z')) then
900 lcase(i:i) = char(ichar(c) + idiff)
905 !-----------------------------------------------------------------------------
906 logical function lcom(ipos,karta)
907 character(len=80) :: karta
908 character :: koment(2) = (/'!','#'/)
913 if (karta(ipos:ipos).eq.koment(i)) lcom=.true.
917 !-----------------------------------------------------------------------------
918 logical function lower_case(ch)
920 lower_case=(ch.ge.'a' .and. ch.le.'z')
922 end function lower_case
923 !-----------------------------------------------------------------------------
924 subroutine mykey(line,keywd,ipos,blankline,errflag)
926 ! This subroutine seeks a non-empty substring keywd in the string LINE.
927 ! The substring begins with the first character different from blank and
928 ! "=" encountered right to the pointer IPOS (inclusively) and terminates
929 ! at the character left to the first blank or "=". When the subroutine is
930 ! exited, the pointer IPOS is moved to the position of the terminator in LINE.
931 ! The logical variable BLANKLINE is set at .TRUE., if LINE(IPOS:) contains
932 ! only separators or the maximum length of the data line (80) has been reached.
933 ! The logical variable ERRFLAG is set at .TRUE. if the string
934 ! consists only from a "=".
935 integer, parameter :: maxlen=80
936 character(len=1) :: empty=' ',equal='=',comma=','
937 character*(*) :: keywd
938 character(len=80) :: line
939 logical :: blankline,errflag !EL,lcom
940 integer :: ipos,istart,iend
942 do while (line(ipos:ipos).eq.empty .and. (ipos.le.maxlen))
945 if (ipos.gt.maxlen .or. lcom(ipos,line) ) then
946 ! At this point the rest of the input line turned out to contain only blanks
947 ! or to be commented out.
953 ! Checks whether the current char is a separator.
954 do while (line(ipos:ipos).ne.empty .and. line(ipos:ipos).ne.equal &
955 .and. line(ipos:ipos).ne.comma .and. ipos.le.maxlen)
959 ! Error flag set to .true., if the length of the keyword was found less than 1.
960 if (iend.lt.istart) then
964 keywd=line(istart:iend)
967 !-----------------------------------------------------------------------------
968 subroutine numstr(inum,numm)
969 character(len=10) :: huj='0123456789'
970 character*(*) :: numm
971 integer :: inumm,inum,inum1,inum2
976 numm(3:3)=huj(inum2+1:inum2+1)
980 numm(2:2)=huj(inum2+1:inum2+1)
984 numm(1:1)=huj(inum2+1:inum2+1)
986 end subroutine numstr
988 !-----------------------------------------------------------------------------
989 function ucase(string)
990 integer :: i, k, idiff
991 character(*) :: string
992 character(len=len(string)) :: ucase
993 character(len=1) :: c
994 character(len=40) :: chtmp
1000 if (string(k+1:) .ne. ' ') then
1004 idiff = ichar('a') - ichar('A')
1008 if (lge(c,'a') .and. lle(c,'z')) then
1009 ucase(i:i) = char(ichar(c) - idiff)
1014 !-----------------------------------------------------------------------------
1016 !-----------------------------------------------------------------------------
1017 subroutine pdbout(etot,tytul,iunit)
1019 use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize
1024 ! implicit real*8 (a-h,o-z)
1025 ! include 'DIMENSIONS'
1026 ! include 'COMMON.CHAIN'
1027 ! include 'COMMON.INTERACT'
1028 ! include 'COMMON.NAMES'
1029 ! include 'COMMON.IOUNITS'
1030 ! include 'COMMON.HEADER'
1031 ! include 'COMMON.SBRIDGE'
1032 ! include 'COMMON.DISTFIT'
1033 ! include 'COMMON.MD'
1034 !el character(len=50) :: tytul
1035 character*(*) :: tytul
1036 character(len=1),dimension(58) :: chainid= (/'A','B','C','D','E','F','G','H','I','J',&
1037 'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',&
1038 'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',&
1039 '0','1','2','3','4','5','6','7','8','9'/)
1041 integer,dimension(maxres) :: ica !(maxres)
1043 integer,dimension(nres) :: ica !(maxres)
1047 integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
1048 real(kind=8) :: etot,xi,yi,zi
1052 if(.not.allocated(molnum)) allocate(molnum(maxres))
1053 ! molnum(:)=mnumi(:)
1054 if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
1055 if(.not.allocated(istype)) allocate(istype(maxres))
1057 if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
1060 write (iunit,'(3a,1pe15.5)') 'REMARK ',tytul,' ENERGY ',etot
1061 !model write (iunit,'(a5,i6)') 'MODEL',1
1063 if (nhfrag.gt.0) then
1065 iti=itype(hfrag(1,j),1)
1066 itj=itype(hfrag(2,j),1)
1068 write (iunit,'(a5,i5,1x,a1,i1,2x,a3,i7,2x,a3,i7,i3,t76,i5)') &
1070 restyp(iti,1),hfrag(1,j)-1,&
1071 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1073 write (iunit,'(a5,i5,1x,a1,i2,1x,a3,i7,2x,a3,i7,i3)') &
1075 restyp(iti,1),hfrag(1,j)-1,&
1076 restyp(itj,1),hfrag(2,j)-1,1,hfrag(2,j)-hfrag(1,j)
1081 if (nbfrag.gt.0) then
1085 iti=itype(bfrag(1,j),1)
1086 itj=itype(bfrag(2,j)-1,1)
1088 write (iunit,'(a5,i5,1x,a1,i1,i3,1x,a3,i6,2x,a3,i6,i3)') &
1090 restyp(iti,1),bfrag(1,j)-1,&
1091 restyp(itj,1),bfrag(2,j)-2,0
1093 if (bfrag(3,j).gt.bfrag(4,j)) then
1095 itk=itype(bfrag(3,j),1)
1096 itl=itype(bfrag(4,j)+1,1)
1098 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)') &
1100 restyp(itl,1),bfrag(4,j),&
1101 restyp(itk,1),bfrag(3,j)-1,-1,&
1102 "N",restyp(itk,1),bfrag(3,j)-1,&
1103 "O",restyp(iti,1),bfrag(1,j)-1
1107 itk=itype(bfrag(3,j),1)
1108 itl=itype(bfrag(4,j)-1,1)
1111 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)') &
1113 restyp(itk,1),bfrag(3,j)-1,&
1114 restyp(itl,1),bfrag(4,j)-2,1,&
1115 "N",restyp(itk,1),bfrag(3,j)-1,&
1116 "O",restyp(iti,1),bfrag(1,j)-1
1128 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1129 'SSBOND',i,'CYS',idssb(i)-nnt+1,&
1130 'CYS',jdssb(i)-nnt+1
1132 write(iunit,'(a6,i4,1x,a3,i7,4x,a3,i7)') &
1133 'SSBOND',i,'CYS',ihpb(i)-nnt+1-nres,&
1134 'CYS',jhpb(i)-nnt+1-nres
1142 write(iout,*) "TUTUT"
1144 ! write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
1145 iti=itype(i,molnum(i))
1147 if (molnum(i+1).eq.0) then
1148 iti1=ntyp1_molec(molnum(i))
1150 iti1=itype(i+1,molnum(i+1))
1152 if ((iti.eq.ntyp1_molec(molnum(i))).and.(iti1.eq.ntyp1_molec(molnum(i)))) cycle
1154 if (iti.eq.ntyp1_molec(molnum(i)).and.(molnum(i+1).eq.5)) cycle
1156 if (iti.eq.ntyp1_molec(molnum(i))) then
1158 if (ichain.gt.58) ichain=1
1160 write (iunit,'(a)') 'TER'
1165 if (molnum(i).eq.1) then
1167 write (iunit,10) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1168 ires,(c(j,i),j=1,3),vtot(i)
1169 elseif(molnum(i).eq.2) then
1170 if (istype(i).eq.0) istype(i)=1
1171 write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
1172 chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
1173 elseif (molnum(i).eq.4) then
1177 ! xi=dmod(xi,boxxsize)
1178 ! if (xi.lt.0.0d0) xi=xi+boxxsize
1179 ! yi=dmod(yi,boxysize)
1180 ! if (yi.lt.0.0d0) yi=yi+boxysize
1181 ! zi=dmod(zi,boxzsize)
1182 ! if (zi.lt.0.0d0) zi=zi+boxzsize
1183 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1184 ires,xi,yi,zi,vtot(i)
1185 elseif (molnum(i).eq.5) then
1189 xi=dmod(xi,boxxsize)
1190 if (xi.lt.0.0d0) xi=xi+boxxsize
1191 yi=dmod(yi,boxysize)
1192 if (yi.lt.0.0d0) yi=yi+boxysize
1193 zi=dmod(zi,boxzsize)
1194 if (zi.lt.0.0d0) zi=zi+boxzsize
1195 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1196 ires,xi,yi,zi,vtot(i)
1199 write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
1200 ires,(c(j,i),j=1,3),vtot(i)
1202 if ((iti.ne.10).and.(molnum(i).ne.5)) then
1204 if (molnum(i).eq.1) then
1205 write (iunit,20) iatom,restyp(iti,1),chainid(ichain),&
1206 ires,(c(j,nres+i),j=1,3),&
1208 else if (molnum(i).eq.2) then
1209 if (istype(i).eq.0) istype(i)=1
1210 write (iunit,50) iatom,sugartyp(istype(i)),restyp(iti,2), &
1211 chainid(ichain),ires,(c(j,nres+i),j=1,3),vtot(i+nres)
1217 write (iunit,'(a)') 'TER'
1219 if (molnum(i).eq.5) cycle
1220 if (itype(i,1).eq.ntyp1) cycle
1221 if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
1222 write (iunit,30) ica(i),ica(i+1)
1223 else if (itype(i,1).ne.10 .and. itype(i+1,1).ne.ntyp1) then
1224 write (iunit,30) ica(i),ica(i+1),ica(i)+1
1225 else if (itype(i,1).ne.10 .and. itype(i+1,1).eq.ntyp1) then
1226 write (iunit,30) ica(i),ica(i)+1
1229 if ((itype(nct,1).ne.10).and.(molnum(i).ne.5)) then
1230 write (iunit,30) ica(nct),ica(nct)+1
1234 write (iunit,30) ica(idssb(i))+1,ica(jdssb(i))+1
1236 write (iunit,30) ica(ihpb(i)-nres)+1,ica(jhpb(i)-nres)+1
1239 write (iunit,'(a6)') 'ENDMDL'
1240 10 FORMAT ('ATOM',I7,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1242 20 FORMAT ('ATOM',I7,' CB ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1243 40 FORMAT ("ATOM",I7," C5' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1244 50 FORMAT ("ATOM",I7," C1' ",1X,2A1,1X,A1,I4,4X,3F8.3,f15.3)
1246 30 FORMAT ('CONECT',8I5)
1247 60 FORMAT ('HETATM',I5,' CA ',A3,1X,A1,I4,4X,3F8.3,f15.3)
1250 end subroutine pdbout
1252 !-----------------------------------------------------------------------------
1253 subroutine MOL2out(etot,tytul)
1254 ! Prints the Cartesian coordinates of the alpha-carbons in the Tripos mol2
1256 use geometry_data, only: c
1259 ! implicit real*8 (a-h,o-z)
1260 ! include 'DIMENSIONS'
1261 ! include 'COMMON.CHAIN'
1262 ! include 'COMMON.INTERACT'
1263 ! include 'COMMON.NAMES'
1264 ! include 'COMMON.IOUNITS'
1265 ! include 'COMMON.HEADER'
1266 ! include 'COMMON.SBRIDGE'
1267 character(len=32) :: tytul,fd
1268 character(len=3) :: zahl
1269 character(len=6) :: res_num,pom !,ucase
1273 real(kind=8) :: etot
1277 #elif (defined CRAY)
1282 write (imol2,'(a)') '#'
1283 write (imol2,'(a)') &
1284 '# Creating user name: unres'
1285 write (imol2,'(2a)') '# Creation time: ',&
1287 write (imol2,'(/a)') '\@<TRIPOS>MOLECULE'
1288 write (imol2,'(a)') tytul
1289 write (imol2,'(5i5)') nct-nnt+1,nct-nnt+nss+1,nct-nnt+nss+1,0,0
1290 write (imol2,'(a)') 'SMALL'
1291 write (imol2,'(a)') 'USER_CHARGES'
1292 write (imol2,'(a)') '\@<TRIPOS>ATOM'
1294 write (zahl,'(i3)') i
1295 pom=ucase(restyp(itype(i,1),1))
1296 res_num = pom(:3)//zahl(2:)
1297 write (imol2,10) i,(c(j,i),j=1,3),i,res_num,0.0
1299 write (imol2,'(a)') '\@<TRIPOS>BOND'
1301 write (imol2,'(i5,2i6,i2)') i-nnt+1,i-nnt+1,i-nnt+2,1
1304 write (imol2,'(i5,2i6,i2)') nct-nnt+i,ihpb(i),jhpb(i),1
1306 write (imol2,'(a)') '\@<TRIPOS>SUBSTRUCTURE'
1308 write (zahl,'(i3)') i
1309 pom = ucase(restyp(itype(i,1),1))
1310 res_num = pom(:3)//zahl(2:)
1311 write (imol2,30) i-nnt+1,res_num,i-nnt+1,0
1313 10 FORMAT (I7,' CA ',3F10.4,' C.3',I8,1X,A,F11.4,' ****')
1314 30 FORMAT (I7,1x,A,I14,' RESIDUE',I13,' **** ****')
1316 end subroutine MOL2out
1317 !-----------------------------------------------------------------------------
1321 use energy_data, only: itype
1323 ! implicit real*8 (a-h,o-z)
1324 ! include 'DIMENSIONS'
1325 ! include 'COMMON.IOUNITS'
1326 ! include 'COMMON.CHAIN'
1327 ! include 'COMMON.VAR'
1328 ! include 'COMMON.LOCAL'
1329 ! include 'COMMON.INTERACT'
1330 ! include 'COMMON.NAMES'
1331 ! include 'COMMON.GEO'
1332 ! include 'COMMON.TORSION'
1336 write (iout,'(/a)') 'Geometry of the virtual chain.'
1337 write (iout,'(7a)') ' Res ',' d',' Theta',&
1338 ' Phi',' Dsc',' Alpha',' Omega'
1341 ! print *,vbld(i),"vbld(i)",i
1342 write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),&
1343 rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
1347 end subroutine intout
1348 !-----------------------------------------------------------------------------
1350 subroutine briefout(it,ener,free)!,plik)
1352 subroutine briefout(it,ener)
1357 ! implicit real*8 (a-h,o-z)
1358 ! include 'DIMENSIONS'
1359 ! include 'COMMON.IOUNITS'
1360 ! include 'COMMON.CHAIN'
1361 ! include 'COMMON.VAR'
1362 ! include 'COMMON.LOCAL'
1363 ! include 'COMMON.INTERACT'
1364 ! include 'COMMON.NAMES'
1365 ! include 'COMMON.GEO'
1366 ! include 'COMMON.SBRIDGE'
1367 ! print '(a,i5)',intname,igeom
1370 real(kind=8) :: ener,free
1371 ! character(len=80) :: plik
1374 #if defined(AIX) || defined(PGI)
1375 open (igeom,file=intname,position='append')
1377 open (igeom,file=intname,access='append')
1386 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1388 WRITE (igeom,180) IT,ENER,free,NSS,(IHPB(I),JHPB(I),I=1,8)
1390 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,NSS)
1392 WRITE (igeom,180) IT,ENER,NSS,(IHPB(I),JHPB(I),I=1,9)
1394 WRITE (igeom,190) (IHPB(I),JHPB(I),I=10,NSS)
1396 ! IF (nvar.gt.nphi) WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1397 WRITE (igeom,200) (RAD2DEG*THETA(I),I=3,NRES)
1398 WRITE (igeom,200) (RAD2DEG*PHI(I),I=4,NRES)
1399 ! if (nvar.gt.nphi+ntheta) then
1400 write (igeom,200) (rad2deg*alph(i),i=2,nres-1)
1401 write (igeom,200) (rad2deg*omeg(i),i=2,nres-1)
1404 180 format (I5,F12.3,I2,9(1X,2I3))
1405 190 format (3X,11(1X,2I3))
1408 end subroutine briefout
1409 !-----------------------------------------------------------------------------
1411 subroutine fdate(fd)
1412 character(len=32) :: fd
1415 end subroutine fdate
1417 !-----------------------------------------------------------------------------
1419 real(kind=8) function gyrate(jcon)
1421 real(kind=8) function gyrate()
1424 use geometry_data, only: c
1427 ! implicit real*8 (a-h,o-z)
1428 ! include 'DIMENSIONS'
1429 ! include 'COMMON.INTERACT'
1430 ! include 'COMMON.CHAIN'
1432 real(kind=8),dimension(3) :: cen
1442 cen(j)=cen(j)+c(j,i)
1446 cen(j)=cen(j)/dble(nct-nnt+1)
1451 rg = rg + (c(j,i)-cen(j))**2
1455 gyrate = dsqrt(rg/dble(nct-nnt+1))
1457 gyrate = sqrt(rg/dble(nct-nnt+1))
1462 !-----------------------------------------------------------------------------
1464 subroutine reads(rekord,lancuch,wartosc,default)
1466 character*(*) :: rekord,lancuch,wartosc,default
1467 character(len=80) :: aux
1468 integer :: lenlan,lenrec,iread,ireade
1472 lenlan=ilen(lancuch)
1474 iread=index(rekord,lancuch(:lenlan)//"=")
1475 ! print *,"rekord",rekord," lancuch",lancuch
1476 ! print *,"iread",iread," lenlan",lenlan," lenrec",lenrec
1477 if (iread.eq.0) then
1481 iread=iread+lenlan+1
1482 ! print *,"iread",iread
1483 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1484 do while (iread.le.lenrec .and. iblnk(rekord(iread:iread)))
1486 ! print *,"|",rekord(iread:iread),"|",iblnk(rekord(iread:iread))
1488 ! print *,"iread",iread
1489 if (iread.gt.lenrec) then
1494 ! print *,"ireade",ireade
1495 do while (ireade.lt.lenrec .and. &
1496 .not.iblnk(rekord(ireade:ireade)))
1499 wartosc=rekord(iread:ireade)
1501 end subroutine reads
1503 !-----------------------------------------------------------------------------
1505 !-----------------------------------------------------------------------------
1506 subroutine permut(isym,nperm,tabperm)
1507 !c integer maxperm,maxsym
1508 !c parameter (maxperm=3628800)
1509 !c parameter (maxsym=10)
1510 ! include "DIMENSIONS"
1511 integer n,a,tabperm,nperm,kkk,i,isym
1514 dimension a(isym),tabperm(50,5040)
1529 !c print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n)
1534 if(nextp(n,a)) go to 10
1539 !-----------------------------------------------------------------------------
1540 logical function nextp(n,a)
1542 integer :: n,i,j,k,t
1544 integer,dimension(n) :: a
1546 10 if(a(i).lt.a(i+1)) go to 20
1563 if(a(j).lt.a(i)) go to 40
1571 !-----------------------------------------------------------------------------
1573 !-----------------------------------------------------------------------------
1574 integer function rescode(iseq,nam,itype,molecule)
1576 ! implicit real*8 (a-h,o-z)
1577 ! include 'DIMENSIONS'
1578 ! include 'COMMON.NAMES'
1579 ! include 'COMMON.IOUNITS'
1580 character(len=3) :: nam !,ucase
1581 integer :: iseq,itype,i
1583 print *,molecule,nam
1584 if (molecule.eq.1) then
1585 if (itype.eq.0) then
1587 do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1588 if (ucase(nam).eq.restyp(i,molecule)) then
1596 do i=-ntyp1_molec(molecule),ntyp1_molec(molecule)
1597 if (nam(1:1).eq.onelet(i)) then
1604 else if (molecule.eq.2) then
1605 do i=1,ntyp1_molec(molecule)
1606 print *,nam(1:1),restyp(i,molecule)(1:1)
1607 print *,nam(2:2),"read",i
1608 if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
1613 else if (molecule.eq.3) then
1614 write(iout,*) "SUGAR not yet implemented"
1616 else if (molecule.eq.4) then
1617 write(iout,*) "Explicit LIPID not yet implemented"
1619 else if (molecule.eq.5) then
1620 do i=-1,ntyp1_molec(molecule)
1621 print *,restyp(i,molecule)
1622 print *,i,restyp(i,molecule)(1:2),nam(1:2)
1623 if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then
1629 write(iout,*) "molecule not defined"
1631 write (iout,10) iseq,nam
1633 10 format ('**** Error - residue',i4,' has an unresolved name ',a3)
1634 end function rescode
1635 integer function sugarcode(sugar,ires)
1638 if (sugar.eq.'D') then
1640 else if (sugar.eq.' ') then
1643 write (iout,*) 'UNKNOWN sugar type for residue',ires,' ',sugar
1647 end function sugarcode
1648 integer function rescode_lip(res,ires)
1649 character(len=3):: res
1652 if (res.eq.'NC3') rescode_lip =4
1653 if (res.eq.'NH3') rescode_lip =2
1654 if (res.eq.'PO4') rescode_lip =3
1655 if (res.eq.'GL1') rescode_lip =12
1656 if (res.eq.'GL2') rescode_lip =12
1657 if (res.eq.'C1A') rescode_lip =18
1658 if (res.eq.'C2A') rescode_lip =18
1659 if (res.eq.'C3A') rescode_lip =18
1660 if (res.eq.'C4A') rescode_lip =18
1661 if (res.eq.'C1B') rescode_lip =18
1662 if (res.eq.'C2B') rescode_lip =18
1663 if (res.eq.'C3B') rescode_lip =18
1664 if (res.eq.'C4B') rescode_lip =18
1665 if (res.eq.'D2A') rescode_lip =16
1666 if (res.eq.'D2B') rescode_lip =16
1668 if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res
1671 !-----------------------------------------------------------------------------
1672 !-----------------------------------------------------------------------------