source: setmvs.f@ 4e219a3

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

Update Lund init and BG/P compatibility.

  • Moved call of init_lundff after initialization of the molecule. The Lund

force field needs the molecule to calculate the hydrophobicity matrix.

  • Added bgp target to Makefile to compile SMMP on a BlueGene/P
  • Split if statement in setmvs to avoid out-of-range warning for array

access.

  • Changed setup of boundaries in Protein.zimmer() to match new version of

numpy.

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

  • Property mode set to 100644
File size: 15.3 KB
Line 
1!**************************************************************
2!
3! This file contains the subroutines: setmvs,fndbrn
4!
5! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
6! Shura Hayryan, Chin-Ku
7! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
8! Jan H. Meinke, Sandipan Mohanty
9!
10! **************************************************************
11
12
13 subroutine setmvs(nml)
14
15! ......................................................................
16! PURPOSE: 1. ORDER variables according to rules:
17! variables with same base: 1st comes TORSION (can be only
18! one with this base, since PHASE a. assumed to be FIXED),
19! after this, for atoms branching from this base:
20! for a b.angle & b.length with common primary moving
21! atom=branch atom - b.angle comes 1st
22!
23! iorvr(i), i=i_fivr_ml,i_lavr_ml -> indices of ordered var.
24!
25! 2. define NON-OVERLAPPING moving sets of atoms in molecule
26! 'nml' related to local variables
27!
28! nmsml(i_ml) - number of moving sets per molecule
29! imsvr1(i_vr),imsvr2() - indices of 1st/last m.s for var. 'i_vr'
30! in 'latms1' & 'latms2'
31! latms1(i_ms),latms2() - range of atoms of i-th m.s
32!
33! 3. define indices of next-following variables for each var.,
34! which complete its physical moving set ('added' variables)
35!
36! nadml(i_ml) - number of 'added' var.s per molecule
37! iadvr1(i_vr),iadvr2() - indices of 1st/last 'added' var. for
38! var. 'i_vr' in 'ladvr'
39! ladvr() - indices of 'added' variables
40!
41! 4. define index of corresponding variable for each atom
42!
43! ! routine must be called successively for molecules 1 -> ntlml
44!
45! CALLS: fndbrn, nursvr
46! ......................................................................
47
48 include 'INCL.H'
49
50 logical bb
51
52 parameter (mxh=10)
53 dimension lvw1h(mxh),lvw2h(mxh),l1h(mxh),l2h(mxh)
54
55
56 ntlvr=nvrml(nml)
57
58 if (nml.eq.1) then
59 imsml1(1)=1
60 nms=0
61 iadml1(1)=1
62 nad=0
63 else
64 imsml1(nml)=imsml1(nml-1)+nmsml(nml-1)
65 nms=imsml1(nml)-1
66 iadml1(nml)=iadml1(nml-1)+nadml(nml-1)
67 nad=iadml1(nml)-1
68 endif
69
70 if (ntlvr.eq.0) then
71 write (*,'(a,i4)')
72 & ' setmvs> No variables defined in molecule #',nml
73 nmsml(nml)=0
74 nadml(nml)=0
75 return
76 endif
77! _________________ Take index of primary atom for each variable
78! (i.e. index of atom moved by variable) to
79! sort variables, handling variables with same base:
80! modify indices to obtain appropriate order
81
82 ifirs=irsml1(nml)
83 ilars=irsml2(nml)
84
85 ifivr=ivrml1(nml)
86 ilavr=ifivr+ntlvr-1
87 ifiat=iatrs1(ifirs)
88 ilaat=iatrs2(ilars)
89
90 do n=ifirs,ilars ! ______________________ Residues
91 ib=ivrrs1(n)
92 do i=ib,ib+nvrrs(n)-1 ! _________________ Variables
93 ia=iatvr(i)
94 io=iowat(ia) ! ('ia' cannot be 1st atom of 'nml')
95 it=ityvr(i)
96 if (it.eq.3) then ! torsion
97 do j=1,nbdat(io)
98 ii=ibdat(j,io)
99 if (iowat(ii).eq.io) ia=min(ia,ii)
100 enddo
101 iadvr1(i)=ia*10
102 elseif (it.eq.2) then ! bond angle
103 iadvr1(i)=ia*10+1
104 elseif (it.eq.1) then ! bond length
105 iadvr1(i)=ia*10+2
106 endif
107 iorvr(i)=i ! (initialize for sorting)
108 enddo ! ... Variables
109 enddo ! ... Residues
110! ___________________________________ Sort variables in ascending order
111! (i.e. from start of molecule/base of branches)
112! array 'iorvr' gives indices of (1st,2nd, ... ,n-th) variables;
113! as can be found in arrays for variables (example: ityvr(iorvr())
114 k=ilavr
115 l=ifivr+ntlvr/2
116 ii=ifivr-1
117 1 if (l.gt.ifivr) then
118 l=l-1
119 io=iorvr(l)
120 n=iadvr1(io)
121 else
122 io=iorvr(k)
123 n=iadvr1(io)
124 iorvr(k)=iorvr(ifivr)
125 k=k-1
126 if (k.eq.ifivr) then
127 iorvr(k)=io
128 goto 2
129 endif
130 endif
131 i=l
132 j=l+l-ii
133 do while (j.le.k)
134! if (j.lt.k.and.iadvr1(iorvr(j)).lt.iadvr1(iorvr(j+1))) j=j+1
135 if (j.lt.k) then
136 if (iadvr1(iorvr(j)).lt.iadvr1(iorvr(j+1))) then
137 j = j + 1
138 end if
139 end if
140 if (n.lt.iadvr1(iorvr(j))) then
141 iorvr(i)=iorvr(j)
142 i=j
143 j=j+j-ii
144 else
145 j=k+1
146 endif
147 enddo
148 iorvr(i)=io
149 goto 1
150! ______________________________ Find non-overlapping ranges of atoms (moving
151! sets) for each variable
152 2 nms=imsml1(nml)-1
153
154 do io=ifivr,ilavr ! _____ Loop over variables in 'ascendent' order
155 iv=iorvr(io)
156 ir=nursvr(iv) ! residue for variable 'iv'
157 ia=iatvr(iv) ! primary mov. atom
158 ib=iowat(ia) ! base
159! __________________________ First, determine complete mov. set for 'iv'
160 it=ityvr(iv)
161 if (it.eq.3) then ! torsion
162 i1=0
163 do i=1,nbdat(ib)
164 j=ibdat(i,ib)
165 if (iowat(j).eq.ib) then ! excl. ring
166 call fndbrn(nml,ir,j,k,irg1,irg2,bb)
167 if (bb) k=ilaat
168 if (i1.ne.0) then ! combine ranges
169 if (j.gt.(i2+1).or.k.lt.(i1-1)) then
170 write (*,'(3a,/,2a,i4,a,i3)')
171 & ' setmvs> Cannot combine disjunct ranges of atom',
172 & ' indices for torsion ',nmvr(iv),' in residue ',
173 & seq(ir),ir,' of molecule # ',nml
174 stop
175 else
176 if (j.lt.i1) i1=j
177 if (k.gt.i2) i2=k
178 endif
179 else
180 i1=j
181 i2=k
182 endif
183 endif
184 enddo
185 elseif (it.eq.2.or.it.eq.1) then ! b. angle, b. length
186 i1=ia
187 call fndbrn(nml,ir,i1,i2,irg1,irg2,bb)
188 if (bb) i2=ilaat
189 endif
190
191 if ((nms+1).gt.mxms) then
192 write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
193 & ': Number of moving sets > ',mxms
194 stop
195 endif
196
197 imsvr1(iv)=nms+1 ! index of 1st
198 imsvr2(iv)=nms+1 ! & last m.s for var. 'iv'
199
200! ______________ Next, exclude overlaps between mov. set for 'iv' and the
201! m.s. for 'previous' variables by reducing/splitting those
202
203 do jo=ifivr,io-1 ! prev. variables ...
204 jv=iorvr(jo)
205
206 j1s=imsvr1(jv) ! index of 1st m.s. for 'jv'
207 jns=imsvr2(jv)-j1s+1 ! # of m.s. for 'jv'
208
209 j=j1s
210 do while (j.lt.(j1s+jns)) ! while there are m.s. for 'jv'
211 j1=latms1(j) ! 1st &
212 j2=latms2(j) ! last atom of m.s. 'j'
213 if (i1.le.j2.and.i2.ge.j1) then ! Overlap
214 ja=0
215 if (i1.gt.j1) then
216 if (i2.gt.j2) goto 6
217 ja=1
218 latms2(j)=i1-1
219 endif
220 if (i2.lt.j2) then
221 if (i1.lt.j1) goto 6
222 if (ja.gt.0) then ! +1 moving set
223 nms=nms+1
224 if (nms.gt.mxms) then
225 write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',
226 & nml,': Number of moving sets > ',mxms
227 stop
228 endif
229 jns=jns+1
230 do k=nms,j+2,-1 ! shift ranges of m.s.
231 latms1(k)=latms1(k-1)
232 latms2(k)=latms2(k-1)
233 enddo
234 do ko=jo+1,io
235 k=iorvr(ko)
236 imsvr1(k)=imsvr1(k)+1
237 imsvr2(k)=imsvr2(k)+1
238 enddo
239 latms2(j+1)=j2
240 endif
241 latms1(j+ja)=i2+1
242 ja=ja+1
243 endif
244 if (ja.eq.0) then ! -1 moving set
245 nms=nms-1
246 jns=jns-1
247 do k=j,nms
248 latms1(k)=latms1(k+1)
249 latms2(k)=latms2(k+1)
250 enddo
251 do ko=jo+1,io
252 k=iorvr(ko)
253 imsvr1(k)=imsvr1(k)-1
254 imsvr2(k)=imsvr2(k)-1
255 enddo
256 else
257 j=j+ja
258 endif
259 else ! No overlap
260 j=j+1
261 endif
262 enddo ! mov. sets for 'jv'
263 imsvr2(jv)=j1s+jns-1
264
265 enddo ! prev. variables
266! _______________________________ Finally, add moving set for 'iv'
267 nms=nms+1
268 latms1(nms)=i1
269 latms2(nms)=i2
270 enddo ! variables
271 nmsml(nml)=nms-imsml1(nml)+1
272! _____________________________ Determine index of moving set for each atom
273 do ia=ifiat,ilaat
274 ixmsat(ia)=0
275 enddo
276 do is=imsml1(nml),nms
277 do ia=latms1(is),latms2(is)
278 ixmsat(ia)=is
279 enddo
280 enddo
281! _____________________________ Determine indices of variables which moving
282! set sets have to be added (=are related) to
283! those of a given variable
284
285 i=iorvr(ifivr) ! initialize index of CURRENT var.
286 ii=imsvr1(i) ! -"- index of its 1st m.s
287
288 do io=ifivr,ilavr-1 ! ________ loop over variables
289
290 ic=i ! save index of CURRENT var.
291 ia=iatvr(i) ! ist primar.mv.atom
292 ib=iowat(ia) ! its base
293 it=ityvr(i) ! its type
294 is=ii ! index of its 1st m.s
295
296 n=nad+1
297 iadvr1(i)=n ! # of its 1st 'added' var.
298
299 i=iorvr(io+1) ! index of next-in-order var.
300 ii=imsvr1(i) ! index of its 1st m.s
301
302 do jo=io+1,ilavr ! ______ over following-in-order var.
303 j=iorvr(jo) ! index of var.
304 ja=iatvr(j) ! its prim.mv.at
305 jb=iowat(ja) ! its base
306
307! _______________ current var. is torsion & shares base with var. 'j'
308 if (it.eq.3.and.jb.eq.ib) then
309 do k=n,nad ! ? has this branch been registered before ?
310 if (iatvr(ladvr(k)).eq.ja) goto 3
311 enddo
312 nad=nad+1
313 if (nad.gt.mxvr) then
314 write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
315 & ': Number of added variables > ',mxvr
316 stop
317 endif
318 ladvr(nad)=j ! save index of 'added' variable
319 endif
320
321 3 if (is.lt.ii) then ! _____ current var. has any m.s:
322 do k=is,ii-1 ! ? base of var. 'j' within m.s ?
323 if (latms1(k).le.jb.and.jb.le.latms2(k)) then
324 do l=n,nad
325 if (iatvr(ladvr(l)).eq.ja) goto 4
326 enddo
327 nad=nad+1
328 if (nad.gt.mxvr) then
329 write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
330 & ': Number of added variables > ',mxvr
331 stop
332 endif
333 ladvr(nad)=j
334 endif
335 4 enddo
336 else ! _____ current var. has no m.s:
337 if (ja.eq.ia) then ! ? share prim.mv.at with var. 'j' ?
338 do k=n,nad
339 if (iatvr(ladvr(k)).eq.ja) goto 5
340 enddo
341 nad=nad+1
342 if (nad.gt.mxvr) then
343 write (*,'(a,i4,a,i5)') ' setmvs> Molecule # ',nml,
344 & ': Number of added variables > ',mxvr
345 stop
346 endif
347 ladvr(nad)=j
348 endif
349 endif
350 5 enddo ! ... following-in-order variables
351 iadvr2(ic)=nad ! last 'added' var. for current var.
352 enddo ! ... variables
353
354 iadvr1(i)=nad+1 ! don't forget last variable
355 iadvr2(i)=nad
356
357 nadml(nml)=nad-iadml1(nml)+1
358! _____________________________________ Summary
359! do io=ilavr,ifivr,-1
360! iv=iorvr(io)
361! ib=iowat(iatvr(iv))
362! i1s=imsvr1(iv)
363! i2s=imsvr2(iv)
364! if (i1s.le.i2s) then
365! do i=i1s,i2s
366! i1=latms1(i)
367! i2=latms2(i)
368! if (i.eq.i1s) then
369! write (*,'(a,i3,7a,i4,3a,i4,a)') 'res # ',nursvr(iv),
370! # ' var: ',nmvr(iv),' base:',nmat(ib),' atoms= ',
371! # nmat(i1),'(',i1,') - ',nmat(i2),'(',i2,')'
372! else
373! write (*,'(39x,2a,i4,3a,i4,a)')
374! # nmat(i1),'(',i1,') - ',nmat(i2),'(',i2,')'
375! endif
376! enddo
377! else
378! write (*,'(a,i3,5a)') 'res # ',nursvr(iv),
379! # ' var: ',nmvr(iv),' base:',nmat(ib),' No atoms'
380! endif
381! i1a=iadvr1(iv)
382! i2a=iadvr2(iv)
383! if (i1a.le.i2a) then
384! write (*,'(a,30(1x,a))') ' Depending variables:',
385! # (nmvr(ladvr(i)),i=i1a,i2a)
386! else
387! write (*,'(a)') ' No dep. variables'
388! endif
389! enddo
390! _____________________________________ Summary - End
391
392 return
393
394 6 write (*,'(a,i4,/,2(a,i5),a)')
395 & ' setmvs> Error in atom numbering of molecule # ',nml,
396 & ': atom ranges for variables # ',iv,' and # ',jv,
397 & ' overlap only PARTLY'
398 stop
399
400 end
401! *******************************************************
402 subroutine fndbrn(nml,nrs,ifirg,ilarg,irg1,irg2,bb)
403
404! .........................................................
405! PURPOSE: determine range [ifirg,ilarg] of atom indices
406! for branch starting from atom 'ifirg' of residue
407! 'nrs' in molecule 'nml'
408! OUTPUT: BB - .t. if 'ifirg' is a backbone atom
409! IRG1 & IRG2 - atom indices of ring-closing bond,
410! if 'ifirg' is INSIDE a ring, but NOT
411! its 1st atom ( in 'multiple' rings
412! only LAST closing bond is given !)
413!
414! CALLS: none
415!
416! .........................................................
417
418 include 'INCL.H'
419
420 logical bb
421 dimension ibd(4)
422
423 ilarg=ifirg
424
425 bb=.false.
426 irg1=0
427
428 ifi=iatrs1(nrs)
429 ila=iatrs2(nrs)
430 ixt=ixatrs(nrs)
431
432 if (ifirg.eq.ifi) then ! = 1st mainchain atom
433 bb=.true.
434 if (nrs.ne.irsml1(nml)) then
435 ilarg=ila
436 else ! 1st residue of 'nml'
437
438 ibd(1)=iowat(ifirg)
439 ibd(2)=ibdat(1,ifirg)
440 ibd(3)=ibdat(2,ifirg)
441 ibd(4)=ibdat(3,ifirg)
442
443 il=0
444 do i=1,nbdat(ifirg)+1
445 ib=ibd(i)
446 if (ib.gt.il.and.iowat(ib).eq.ifirg) il=ib
447 enddo
448 if (il.gt.0) ilarg=il-1
449 endif
450 else
451 if (ifirg.eq.ixt) bb=.true.
452 do i=1,nbdat(ifirg) ! ______________ check bonds
453 ib=ibdat(i,ifirg)
454 if (iowat(ib).eq.ifirg) then ! branch
455 do j=ib,ila
456 if (j.gt.ib.and.iowat(j).lt.ib) goto 1
457 if (j.eq.ixt) bb=.true.
458 do k=1,nbdat(j)
459 jb=ibdat(k,j)
460 if (jb.lt.ifirg) then ! ring
461 irg1=j
462 irg2=jb
463 endif
464 enddo
465 ilarg=j
466 enddo ! ... branch atoms
467 elseif (ib.lt.ifirg) then ! ring
468 irg1=ifirg
469 irg2=ib
470 endif
471 1 enddo ! ... bonds
472 endif
473
474 return
475 end
476
Note: See TracBrowser for help on using the repository browser.