source: addend.f@ 32289cd

Last change on this file since 32289cd was cb47b9c, checked in by baerbaer <baerbaer@…>, 14 years ago

Explicitly declare variables.

All variables should be declared so that we can remove the implicit statements
from the beginning of the INCL.H file.

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

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