source: addend.f@ bd2278d

Last change on this file since bd2278d was bd2278d, checked in by baerbaer <baerbaer@…>, 16 years ago

Reformatting comments and continuation marks.

Fortran 90 and higher use ! to mark comments no matter where they are in the
code. The only valid continuation marker is &.
I also added the SMMP.kdevelop.filelist to the repository to make it easier
to use kdevelop.

git-svn-id: svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@12 26dc1dd8-5c4e-0410-9ffe-d298b4865968

  • Property mode set to 100644
File size: 21.3 KB
Line 
1! **************************************************************
2!
3!
4! This file contains the subroutines: addend, redchg, rplgrp
5!
6! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
7! Shura Hayryan, Chin-Ku
8! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
9! Jan H. Meinke, Sandipan Mohanty
10!
11! $Id: addend.f 334 2007-08-07 09:23:59Z meinke $
12! **************************************************************
13 subroutine addend(nml,grpn,grpc)
14
15! ..............................................................
16! PURPOSE: modify terminal residues to complete bonding scheme
17! with residue 'grpn' at N- and residue 'grpc' at C-terminus
18! ! need initial co-ordinates for residues to modify
19! ! for N-terminus: may add only simple groups
20!
21! CALLS: rplgrp,tolost,redchg
22! ..............................................................
23
24 include 'INCL.H'
25
26 character grpn*4,grpc*4
27
28 character res*4,rpat*4,sbrs*4,grn*4,grc*4
29
30
31 grn = grpn
32 call tolost(grn)
33 grc = grpc
34 call tolost(grc)
35
36 if (grn(:3).eq.'ace'.or.grc(:3).eq.'ace'
37 &.or.grn(:3).eq.'nme'.or.grc(:3).eq.'nme') then
38
39 write(*,'(2a)') ' addend> N-Acetyl (ace) or N-Methylamide (nme)'
40 & ,' should be put in SEQUENCE file, not added as end groups'
41
42 stop
43 endif
44
45! __________________________________________ N-terminus
46 ifirs=irsml1(nml)
47 rpat='n '
48 res=seq(ifirs)
49 call tolost(res)
50 if (res(:3).ne.'ace') then
51 if (grn(:3).eq.'nh2') then
52 if (res(:3).eq.'hyp'.and..not.flex) then
53 if(sh2) then
54 sbrs='nh2+'
55 else
56 write (*,'(2a)') ' addend> ',
57 & ' No N-terminal Hyp possible with ECEPP/3 dataset'
58 stop
59 endif
60
61 elseif (res(:3).eq.'pro') then
62 sbrs='nh1 '
63 elseif (res(:3).eq.'gly'.and..not.flex) then
64 if (sh2) then
65 sbrs='nh2 '
66 else
67 sbrs='nh2g'
68 endif
69 else
70 sbrs='nh2 '
71 endif
72
73 elseif (grn.eq.'nh3+') then
74 if (res(:3).eq.'pro'.and..not.flex) then
75 sbrs='nh2+'
76 elseif (res(:3).eq.'gly'.and..not.flex) then
77 sbrs='nh3g'
78 else
79 sbrs='nh3+'
80 endif
81
82 else
83
84 write(*,'(2a)') ' addend> Can add only ',
85 & 'nh2 or nh3+ to N-terminus'
86 stop
87
88 endif
89
90 call rplgrp(nml,ifirs,rpat,sbrs)
91 if (flex) call redchg(nml,ifirs,rpat,sbrs) ! Flex dataset
92
93 else ! ace
94
95 write(*,'(2a)') ' addend> Acetyl group',
96 & ' at N-terminus not modified'
97 endif
98
99! __________________________________________ C-terminus
100 ilars=irsml2(nml)
101 rpat='c '
102 res=seq(ilars)
103 call tolost(res)
104
105 if (res(:3).ne.'nme') then
106
107 if (grc.eq.'cooh') then
108
109 if (res(:3).eq.'gly'.and..not.sh2) then
110 sbrs='gooh'
111 else
112 sbrs='cooh'
113 endif
114
115 elseif (grc.eq.'coo-') then
116
117 if (res(:3).eq.'gly'.and..not.sh2) then
118 sbrs='goo-'
119 else
120 sbrs='coo-'
121 endif
122
123 else
124
125 write(*,'(2a)') ' addend> Can add only ',
126 & 'cooh or coo- to C-terminus'
127 stop
128
129 endif
130
131 call rplgrp(nml,ilars,rpat,sbrs)
132 if (flex) call redchg(nml,ilars,rpat,sbrs) ! Flex dataset
133
134 else ! N'-methylamide
135
136 write(*,'(2a)') ' addend> N-Methylamide',
137 & ' at C-terminus not modified'
138
139 endif
140
141! ----------------------------- net charge of molecule
142 cg = 0.d0
143 do i=iatrs1(irsml1(nml)),iatrs2(irsml2(nml))
144 cg = cg + cgat(i)
145 enddo
146 if (abs(cg).gt.1.d-5) write(*,'(a,i2,a,f7.3,/)')
147 & ' addend> Net charge of molecule #'
148 & ,nml,': ',cg
149
150 return
151 end
152! ****************************************
153 subroutine rplgrp(nml,nrs,rpat,sbrs)
154
155! ...............................................................
156! PURPOSE: replace atom(s) rooted at atom 'rpat' in residue
157! 'nrs' of molecule 'nml' by atom(s) rooted at
158! 'rpat' of residue 'sbrs' (same name of root
159! atom 'rpat' maintains bonding geometry for
160! preceeding atoms in 'nrs')
161!
162! is NOT performed if 'rpat' is within mainchain,
163! except it is first/last mainchain atom of 'nml'
164!
165! CALLS: dihedr,iopfil,iendst,eyring,fndbrn,redres,setsys,valang
166! ...............................................................
167
168 include 'INCL.H'
169
170 character rpat*4,sbrs*4
171 logical ntbb,bb
172
173 dimension ibd(mxbd+1),iybd(mxbd+1)
174
175
176 ifirs=irsml1(nml)
177 ilars=irsml2(nml)
178
179 nxt=ixatrs(nrs)
180 nfi=iatrs1(nrs)
181 nla=iatrs2(nrs)
182! __________________________ indices of atoms to be replaced
183 do i=nfi,nla
184 if (rpat.eq.nmat(i)) then
185 nfirp=i
186 goto 1
187 endif
188 enddo
189 write (*,'(4a,i4,a,i4)') ' rplgrp> cannot find atom >',rpat,
190 &'< to be replaced in residue ',seq(nrs),nrs,' of molecule ',nml
191 stop
192
193 1 call fndbrn(nml,nrs,nfirp,nlarp,irng1,irng2,bb)
194 if (irng1.gt.0) goto 10
195 ntbb=.false.
196
197 if (bb) then ! ..... backbone
198
199 if (nfirp.eq.nxt.and.nrs.eq.ilars) goto 2 ! last mainchain atom
200
201 if (nfirp.eq.nfi.and.nrs.eq.ifirs) then ! 1st MAINCHAIN ATOM !
202
203 ntbb=.true.
204
205 ibd(1)=iowat(nfirp)
206 iybd(1)=iyowat(nfirp)
207
208 do i=1,mxbd
209 ibd(i+1)=ibdat(i,nfirp)
210 iybd(i+1)=iybdat(i,nfirp)
211 enddo
212
213 ibdrg=0 ! __________________ check ring
214 iybdrg=1
215
216 do i=1,nbdat(nfirp)+1
217 if (iowat(ibd(i)).ne.nfirp) then
218 if (ibdrg.ne.0) then
219 write (*,'(2a,i3)')
220 & ' rplgrp> Can handle only simple ring at 1st',
221 & ' atom of molecule #',nml
222 stop
223 endif
224 ibdrg=ibd(i)
225 iybdrg=iybd(i)
226 endif
227 enddo
228 nxtbb1=nlarp+1 ! _________________ next backbone atoms
229 isgbb1=iyowat(nxtbb1)
230
231 if (nxtbb1.eq.nxt) then
232 if (nrs.lt.ilars) then
233 nxtbb2=iatrs1(nrs+1)
234 goto 3
235 endif
236 else
237 do i=nxt,nxtbb1+1,-1
238 if (iowat(i).eq.nxtbb1) then
239 nxtbb2=i
240 goto 3
241 endif
242 enddo
243 endif
244 goto 11
245 else
246 write (*,'(4a,i4,a,i4)')
247 & ' rplgrp> Cannot replace BACKBONE atom ',rpat,
248 & ' of residue ',seq(nrs),nrs,' in molecule #',nml
249 stop
250 endif
251
252 endif ! N-terminus
253! _________________________________ previous atoms
254 2 if (nfirp.eq.nfi.and.nrs.eq.ifirs) goto 11
255 nxtbb1=iowat(nfirp)
256 if (nxtbb1.eq.nfi.and.nrs.eq.ifirs) goto 11
257 nxtbb2=iowat(nxtbb1)
258! _______________________________ get data for substituent atoms
259 3 if (iopfil(lunlib,reslib,'old','formatted').le.izero) then
260 write (*,'(a,/,a,i3,2a)')
261 & ' rplgrp> ERROR opening library of residues:',
262 & ' LUN=',lunlib,' FILE=',reslib(1:iendst(reslib))
263 stop
264 endif
265 call redres(sbrs,natsb,nxtsb,nvrsb)
266 close (lunlib)
267! __________________________ indices of substituent atoms
268 do i=1,natsb
269 if (rpat.eq.nmath(i)) then
270 nfisb=i
271 goto 4
272 endif
273 enddo
274 write (*,'(4a)') ' rplgrp> Cannot find atom >',rpat,
275 &'< in substituent residue ',sbrs
276 stop
277
278 4 nlasb=nfisb
279 do i=1,nbdath(nfisb)
280 ib=ibdath(i,nfisb)
281 if (iowath(ib).lt.nfisb) goto 10
282 do j=ib,natsb
283 if (j.gt.ib.and.iowath(j).lt.ib) goto 5
284 do k=1,nbdath(j)
285 if (ibdath(k,j).lt.nfisb) goto 10
286 enddo
287 nlasb=j
288 enddo ! ... branch atoms
289 5 enddo ! ... branches
290! _________________________________________________ local axes at 'nfirp'
291 call setsys(nxtbb1,nfirp,nxtbb2,x1,x2,x3,y1,y2,y3,z1,z2,z3)
292
293 xtoat(nfirp)=x1
294 ytoat(nfirp)=x2
295 ztoat(nfirp)=x3
296 xbaat(nfirp)=z1
297 ybaat(nfirp)=z2
298 zbaat(nfirp)=z3
299
300! _____________________ add virtual atoms
301 if (ntbb) then
302
303 ct=cstoat(nxtbb2) ! t.angle_(+2)
304 st=sntoat(nxtbb2)
305 ca=csbaat(nxtbb1) ! b.angle_(+1)
306 sa=snbaat(nxtbb1)
307
308! ------------------- Eyring
309 h2=-sa*ct
310 h3=-sa*st
311 x1=-ca*x1+h2*y1+h3*z1
312 x2=-ca*x2+h2*y2+h3*z2
313 x3=-ca*x3+h2*y3+h3*z3
314 dx=one/sqrt(x1*x1+x2*x2+x3*x3)
315 x1=x1*dx
316 x2=x2*dx
317 x3=x3*dx
318
319 xat(izero)=xat(nfirp)+x1
320 yat(izero)=yat(nfirp)+x2
321 zat(izero)=zat(nfirp)+x3
322 z1=-st*y1+ct*z1
323 z2=-st*y2+ct*z2
324 z3=-st*y3+ct*z3
325 dz=one/sqrt(z1*z1+z2*z2+z3*z3)
326 z1=z1*dz
327 z2=z2*dz
328 z3=z3*dz
329
330 ct=cstoat(nxtbb1) ! t.angle_(+1)
331 st=sntoat(nxtbb1)
332
333! -------------------- Eyring with b.angle = 90 deg.
334 xat(-ione)=xat(izero)-ct*(z2*x3-z3*x2)-st*z1
335 yat(-ione)=yat(izero)-ct*(z3*x1-z1*x3)-st*z2
336 zat(-ione)=zat(izero)-ct*(z1*x2-z2*x1)-st*z3
337
338 endif
339! _____________________________________________ Shift atom data
340 nrp=nlarp-nfirp
341 nsb=nlasb-nfisb
342 if (nrp.ne.nsb) then
343 nsh=nsb-nrp
344 do i=nfi,nfirp-1 ! bonds to atoms after repl. group
345 do j=1,mxbd
346 ib=ibdat(j,i)
347 if (ib.gt.nlarp.and.ib.le.nla) ibdat(j,i)=ib+nsh
348 enddo
349 enddo
350 ilaat=iatrs2(irsml2(ntlml))
351 if (nrp.gt.nsb) then ! less atoms
352 i1=nlarp+1
353 i2=ilaat
354 i3=1
355 else ! more atoms
356 if ((ilaat+nsh).gt.mxat) then
357 write (*,'(a,i5)') ' rplgrp> number of atoms > ',mxat
358 stop
359 endif
360
361 i1=ilaat
362 i2=nlarp+1
363 i3=-1
364 endif
365
366 do i=i1,i2,i3
367 ii=i+nsh
368 nbdat(ii)=nbdat(i)
369
370 ibd(1)=iowat(i)
371 iybd(1)=iyowat(i)
372 do j = 1,mxbd
373 ibd(j+1)=ibdat(j,i)
374 iybd(j+1)=iybdat(j,i)
375 enddo
376
377 do j=1,mxbd+1
378 if (ibd(j).gt.nfirp) ibd(j)=ibd(j)+nsh
379 enddo
380
381 iowat(ii)=ibd(1)
382 iyowat(ii)=iybd(1)
383 do j = 1,mxbd
384 ibdat(j,ii)=ibd(j+1)
385 iybdat(j,ii)=iybd(j+1)
386 enddo
387
388 ityat(ii)=ityat(i)
389 nmat(ii)=nmat(i)
390 cgat(ii)=cgat(i)
391 blat(ii)=blat(i)
392 baat(ii)=baat(i)
393 snbaat(ii)=snbaat(i)
394 csbaat(ii)=csbaat(i)
395 toat(ii)=toat(i)
396 sntoat(ii)=sntoat(i)
397 cstoat(ii)=cstoat(i)
398 xat(ii)=xat(i)
399 yat(ii)=yat(i)
400 zat(ii)=zat(i)
401 xtoat(ii)=xtoat(i)
402 ytoat(ii)=ytoat(i)
403 ztoat(ii)=ztoat(i)
404 xbaat(ii)=xbaat(i)
405 ybaat(ii)=ybaat(i)
406 zbaat(ii)=zbaat(i)
407
408 enddo
409! ____________________________________________ Shift residue data
410 do i=nrs+1,irsml2(ntlml)
411 iatrs1(i)=iatrs1(i)+nsh
412 iatrs2(i)=iatrs2(i)+nsh
413 ixatrs(i)=ixatrs(i)+nsh
414 enddo
415 iatrs2(nrs)=nla+nsh
416 if (nxt.gt.nlarp) ixatrs(nrs)=nxt+nsh
417 else
418 nsh=0
419 endif
420! _________________________________________ Correct data of 'nfirp'
421 ish=nfirp-nfisb
422 ityat(nfirp)=ityath(nfisb)
423
424 if(.not.sh2) cgat(nfirp)=cgath(nfisb) ! NOT for 'ECEPP/2'
425
426 nb=nbdath(nfisb) ! _______________ Bonds
427
428 do i=1,mxbd
429 ib=ibdath(i,nfisb)
430 if (ib.ne.0) then
431 ibd(i)=ib+ish
432 iybd(i)=iybdath(i,nfisb)
433 else
434 ibd(i)=0
435 iybd(i)=1
436 endif
437
438 enddo
439
440 if (ntbb) then
441
442 ibd(mxbd+1)=0
443 iybd(mxbd+1)=1
444
445 ibd(nb+1)=nxtbb1+nsh ! bond to next backbone atom
446 iybd(nb+1)=isgbb1
447
448 if (ibdrg.ne.0) then
449 nb=nb+1
450 if (nb.gt.mxbd) then
451 write (*,'(6a,/,2a,3(i4,a))')
452 & ' rplgrp> Cannot add atoms following ',rpat,
453 & ' from group ',sbrs,' to atom ',rpat,
454 & ' of residue ',seq(nrs),nrs,' in molecule #',nml,
455 & ' because need >',(mxbd+1),' bonds'
456 stop
457 endif
458 ibd(nb+1)=ibdrg+nsh
459 iybd(nb+1)=iybdrg
460 endif
461 iowat(nfirp)=ibd(1)
462 iyowat(nfirp)=iybd(1)
463 do i=1,mxbd
464 ibdat(i,nfirp)=ibd(i+1)
465 iybdat(i,nfirp)=iybd(i+1)
466 enddo
467 else ! not 'ntbb'
468
469 do i=1,mxbd
470 ibdat(i,nfirp)=ibd(i)
471 iybdat(i,nfirp)=iybd(i)
472 enddo
473
474 endif
475 nbdat(nfirp)=nb
476! _________________________________________ Add data for substituent
477 ii=nfirp
478 do i=nfisb+1,nlasb
479 ii=ii+1
480 nbdat(ii)=nbdath(i)
481
482 iow=iowath(i)+ish
483 iowat(ii)=iow
484 iyowat(ii)=iyowath(i)
485
486 do j=1,mxbd
487 ib=ibdath(j,i)
488 if (ib.ge.nfisb) then
489 ibdat(j,ii)=ib+ish
490 else
491 ibdat(j,ii)=ib
492 endif
493 iybdat(j,ii)=iybdath(j,i)
494 enddo
495
496 ityat(ii)=ityath(i)
497 nmat(ii)=nmath(i)
498 cgat(ii)=cgath(i)
499 blat(ii)=blath(i)
500 ba=baath(i)
501 baat(ii)=ba
502 csbaat(ii)=cos(ba)
503 snbaat(ii)=sin(ba)
504 to=toath(i)
505 toat(ii)=to
506 cstoat(ii)=cos(to)
507 sntoat(ii)=sin(to)
508 call eyring(ii,iow)
509 if (ntbb) then ! reset some internal coordinates
510 if (iow.eq.nfirp) then
511 ba=valang(izero,nfirp,ii)
512 baat(ii)=ba
513 csbaat(ii)=cos(ba)
514 snbaat(ii)=sin(ba)
515 to=dihedr(-ione,izero,nfirp,ii)
516 toat(ii)=to
517 cstoat(ii)=cos(to)
518 sntoat(ii)=sin(to)
519 elseif (iowat(iow).eq.nfirp) then
520 to=dihedr(izero,nfirp,iow,ii)
521 toat(ii)=to
522 cstoat(ii)=cos(to)
523 sntoat(ii)=sin(to)
524 endif
525
526 endif ! ntbb
527
528 enddo ! substituent atoms
529! ___________________________________________________ Take care of Variables
530! (assume variables of replaced group/substituent to be stored CONSECUTIVELY)
531
532 ilavr=ivrml1(ntlml)+nvrml(ntlml)-1
533 ifivr=ivrrs1(nrs) ! variables to be replaced (#=ivrrp,last=lvrrp)
534 ivrrp=0
535 do i=ifivr,ifivr+nvrrs(nrs)-1
536 iat=iatvr(i)
537 if (nfirp.lt.iat.and.iat.le.nlarp) then
538 ivrrp=ivrrp+1
539 lvrrp=i
540 endif
541 enddo
542 if (ivrrp.eq.0) then ! No variables to replace
543 do i=ifivr,ilavr
544 if (iatvr(i).gt.nlarp) then
545 lvrrp=i-1
546 goto 6
547 endif
548 enddo
549 lvrrp=ilavr
550 endif
551 6 ivrsb=0 ! variables from substituent (#=ivrsb)
552 do i=1,nvrsb
553 iat=iatvrh(i)
554 if (nfisb.lt.iat.and.iat.le.nlasb) then
555 ity=ityvrh(i)
556 if (ntbb.and.iowath(iat).eq.nfisb.and.ity.gt.2) goto 7
557 ivrsb=ivrsb+1
558 nmvrh(ivrsb)=nmvrh(i)
559 ityvrh(ivrsb)=ityvrh(i)
560 iclvrh(ivrsb)=iclvrh(i)
561 iatvrh(ivrsb)=iat
562 endif
563 7 enddo
564 if (nsh.ne.0) then ! if # of atoms changed
565 do i=lvrrp+1,ilavr
566 iatvr(i)=iatvr(i)+nsh
567 enddo
568 endif
569
570 if (ivrrp.ne.ivrsb) then ! shift data for variables
571 jsh=ivrsb-ivrrp
572 if (ivrrp.gt.ivrsb) then
573 i1=lvrrp+1
574 i2=ilavr
575 i3=1
576 else
577 if ((ilavr+jsh).gt.mxvr) then
578 write (*,'(a,i5)') ' rplgrp> number of variables > ',mxvr
579 stop
580 endif
581 i1=ilavr
582 i2=lvrrp+1
583 i3=-1
584 endif
585 do i=i1,i2,i3
586 ii=i+jsh
587 ityvr(ii)=ityvr(i)
588 iclvr(ii)=iclvr(i)
589 nmvr(ii)=nmvr(i)
590 iatvr(ii)=iatvr(i)
591 enddo
592
593 do i=nrs+1,irsml2(ntlml)
594 ivrrs1(i)=ivrrs1(i)+jsh
595 enddo
596 nvrrs(nrs)=nvrrs(nrs)+jsh
597 nvrml(nml)=nvrml(nml)+jsh
598 do i=nml+1,ntlml
599 ivrml1(i)=ivrml1(i)+jsh
600 enddo
601 endif
602 ii=lvrrp-ivrrp ! Add variables for substitutent
603 do i=1,ivrsb
604 ii=ii+1
605 nmvr(ii)=nmvrh(i)
606 ityvr(ii)=ityvrh(i)
607 iclvr(ii)=iclvrh(i)
608 iatvr(ii)=iatvrh(i)+ish
609 enddo
610
611 return
612! __________________________________________ Errors
613 10 write (*,'(3a,/,2a,i4,a,i4,/,2a)')
614 & ' rplgrp> Cannot replace atom(s) following ',rpat,
615 & ' from INSIDE a ring',' in residue: ',seq(nrs),nrs,
616 & ' in molecule #',nml,' or in substitute: ',sbrs
617 stop
618 11 write (*,'(4a,i4,a,i4,/,a)')
619 & ' rplgrp> Cannot replace atom(s) following ',rpat,
620 & ' of residue ',seq(nrs),nrs,' in molecule #',nml,
621 & ' since necessary 2 previous atoms are not available'
622 stop
623
624 end
625! ****************************************
626 subroutine redchg(nml,nrs,rpat,sbrs)
627
628! .........................................................
629! PURPOSE: read and place atomic point charges from residue
630! 'sbrs' to residue 'nrs' of molecule 'nml'
631! from library 'chglib' with LUN=lunchg, if ilib=1
632! 'reslib' with LUN=lunlib, if ilib=2
633!
634! CALLS: iopfil,iendst,tolost
635! ........................................................
636
637 include 'INCL.H'
638
639 character rpat*4,sbrs*4,atnm*4,res*4,cgty*5,line*100
640
641 ifirs=irsml1(nml)
642 ilars=irsml2(nml)
643 res=seq(nrs)
644 call tolost(res)
645
646 if (nrs.eq.ifirs.or.nrs.eq.ilars) then
647
648 if (rpat.eq.'n '.or.rpat.eq.'c ') then
649 if (nrs.eq.ifirs.and.nrs.eq.ilars) goto 10 ! Dont have this yet
650
651 if (rpat.eq.'n '.and.nrs.eq.ifirs) then
652 if (res(1:3).eq.'pro'.or.res(1:3).eq.'hyp') then
653 if (sbrs.eq.'nh1 ') cgty='n'//res
654 if (sbrs.eq.'nh2+') cgty='+'//res
655 else
656 if (sbrs.eq.'nh2 ') cgty='n'//res
657 if (sbrs.eq.'nh3+') cgty='+'//res
658 endif
659 elseif (rpat.eq.'c '.and.nrs.eq.ilars) then
660 if (sbrs.eq.'cooh') cgty='c'//res
661 if (sbrs.eq.'coo-') cgty='-'//res
662 else
663 goto 10
664 endif
665 else
666 write (*,*) ' redchg> dont know which end goup is present'
667 stop
668 endif
669 ilib=1
670 else
671 cgty(1:)=sbrs
672 ilib=2
673 endif
674
675 if (ilib.eq.1) then
676
677 if (iopfil(lunchg,chgfil,'old','formatted').le.izero) then
678 write (*,'(a,/,a,i3,2a)')
679 & ' redchg> ERROR opening library of charges:',
680 & ' LUN=',lunchg,' FILE=',chgfil(1:iendst(chgfil))
681 stop
682 endif
683
684 1 line=' '
685 read (lunchg,'(a)',end=3) line
686 l=iendst(line)
687
688 if (l.ge.10.and.line(1:1).eq.'#'.and.line(2:6).eq.cgty) then
689 read (line(7:10),'(i4)') nchg
690 if (nchg.le.mxath) then
691 read (lunchg,'(10(2x,a4,1x))') (nmath(i),i=1,nchg)
692 read (lunchg,'(10(f6.4,1x))') (cgath(i),i=1,nchg)
693 close(lunchg)
694 do i=iatrs1(nrs),iatrs2(nrs)
695 atnm=nmat(i)
696 do j=1,nchg
697 if (nmath(j).eq.atnm) then
698 cgat(i)=cgath(j)
699 goto 2
700 endif
701 enddo
702 write (*,'(6a)') ' redchg> Cannot find atom: ',atnm,
703 & ' for entry: ',cgty,' in library: ',
704 & chgfil(1:iendst(chgfil))
705 stop
706 2 enddo
707 return
708 else
709 write (*,'(4a)')
710 & ' redchg> must increase MXATH to read data for entry: ',
711 & cgty,' in library: ',chgfil(1:iendst(chgfil))
712 close(lunchg)
713 stop
714 endif
715 endif
716 goto 1
717 3 write (*,'(4a)')
718 & ' redchg> Cannot find entry: ',cgty,' in library: ',
719 & chgfil(1:iendst(chgfil))
720 close(lunchg)
721 stop
722
723 elseif (ilib.eq.2) then
724
725 if (iopfil(lunlib,reslib,'old','formatted').le.izero) then
726 write (*,'(a,/,a,i3,2a)')
727 & ' redchg> ERROR opening library of residues:',
728 & ' LUN=',lunlib,' FILE=',reslib(1:iendst(reslib))
729 stop
730 endif
731
732 4 line=' '
733 read (lunlib,'(a)',end=6) line
734 l=iendst(line)
735
736 if (l.ge.9.and.line(1:1).eq.'#'.and.line(2:5).eq.cgty(1:4)) then
737 read (line(6:9),'(i4)') nchg
738 if (nchg.le.mxath) then
739 read (lunlib,'(a4,42x,d7.4)') (nmath(i),cgath(i),i=1,nchg)
740 close(lunlib)
741 do i=iatrs1(nrs),iatrs2(nrs)
742 atnm=nmat(i)
743 do j=1,nchg
744 if (nmath(j).eq.atnm) then
745 cgat(i)=cgath(j)
746 goto 5
747 endif
748 enddo
749 write (*,'(6a)') ' redchg> Cannot find atom: ',atnm,
750 & ' for entry: ',cgty,' in library: ',
751 & reslib(1:iendst(reslib))
752 stop
753 5 enddo
754 return
755 else
756 write (*,'(4a)')
757 & ' redchg> must increase MXATH to read data for entry: ',
758 & cgty,' in library: ',reslib(1:iendst(reslib))
759 close(lunchg)
760 stop
761 endif
762 endif
763 goto 4
764 6 write (*,'(4a)')
765 & ' redchg> Cannot find entry: ',cgty,' in library: ',
766 & reslib(1:iendst(reslib))
767 close(lunchg)
768 stop
769
770 endif
771
772 10 write (*,'(4a)')
773 & ' redchg> Do not have charges for N/C-terminal residue ',
774 & res,' modified with group :',sbrs
775 stop
776
777 end
778
Note: See TracBrowser for help on using the repository browser.