1      logical function argos_prep_mknoe(lfnout,lfntop,
2     + filtop,lfncmd,filcmd,
3     + lfnnoe,filnoe,mnoe,lfnpmf,filpmf,slvnam)
4c
5c $Id$
6c
7      implicit none
8c
9#include "mafdecls.fh"
10#include "util.fh"
11c
12      logical argos_prep_topsiz,argos_prep_toprd,
13     + argos_prep_disres,argos_prep_select,argos_prep_equiv
14      external argos_prep_topsiz,argos_prep_toprd,
15     + argos_prep_disres,argos_prep_select,argos_prep_equiv
16c
17      integer lfnout,lfntop,lfncmd,lfnnoe,lfnpmf
18      character*255 filtop,filcmd,filnoe,filpmf
19      character*255 card,card2
20c
21      integer nsa,nwa,noe
22      integer mat,msa,mwa,mnoe,msm,nsm,msb,nsb
23      integer nnoe,nqu,mseq,nseq
24      character*3 slvnam
25c
26      integer l_cwa,i_cwa,l_mas,i_mas,l_num,i_num,l_sat,i_sat
27      integer l_csa,i_csa,l_sgm,i_sgm,l_sml,i_sml,l_qsa,i_qsa
28      integer l_inoe,i_inoe,l_dnoe,i_dnoe,l_qu,i_qu
29      integer l_qwa,i_qwa,l_wat,i_wat,l_sfr,i_sfr
30      integer l_isb,i_isb,l_csb,i_csb,l_eq,i_eq
31      integer l_temp,i_temp,l_tempj,i_tempj
32      integer l_lseq,i_lseq,i_ihop,l_ihop,l_istat,i_istat
33c
34      integer maxgrp,maxatm,maxpmf
35      parameter(maxgrp=500)
36      parameter(maxatm=2500)
37      parameter(maxpmf=500)
38c
39      integer igroup(maxgrp,maxatm),lgroup(maxgrp),mgroup(maxgrp)
40      integer numpmf,ipmf(maxpmf,5),iopt(maxpmf,3)
41      real*8 rpmf(maxpmf,4)
42      integer npmfr,ipmfr,nsarep
43c
44      integer i
45c
46      npmfr=0
47      ipmfr=0
48c
49      do 2 i=1,maxgrp
50      lgroup(i)=0
51      mgroup(i)=0
52    2 continue
53      numpmf=0
54c
55      if(util_print('where',print_debug)) then
56      write(lfnout,1000)
57 1000 format(' NOE FILE GENERATION')
58      endif
59c
60      if(.not.argos_prep_topsiz(lfntop,filtop,lfnout,
61     + mat,msa,mwa,msb,nqu,nseq))
62     + call md_abort('argos_prep_topsiz failed',9999)
63c
64      if(util_print('where',print_high)) then
65      write(lfnout,1001) filtop(1:index(filtop,' ')-1),mat,msa,mwa,msb
66 1001 format(' Topology',t40,a,//,
67     + ' Number of atom types',t40,i8,/,
68     + ' Number of solute atoms',t40,i8,/,
69     + ' Number of solvent atoms',t40,i8,/,
70     + ' Number of solute bonds',t40,i8,/)
71      endif
72      nwa=0
73      nsa=0
74      nsb=0
75      nnoe=0
76      mseq=nseq
77c
78c     allocate memory
79c     ---------------
80c
81c     character*16 cwa(mwa)  : solvent atom names
82c     character*16 csa(msa)  : solute atom names
83c     integer isar(msa)      : solute atom types
84c     integer isgm(msa)      : solute segment numbers
85c     integer isfnd(msa)     : solute atom found flags
86c     real*8 qsa(msa)        : solute atom charges
87c     real*8 xs(3,msa)       : solute atom coordinates
88c
89      if(.not.ma_push_get(mt_dbl,nqu,'qu',l_qu,i_qu))
90     + call md_abort('Memory allocation failed for qu',9999)
91      if(.not.ma_push_get(mt_int,2*mnoe,'inoe',l_inoe,i_inoe))
92     + call md_abort('Memory allocation failed for inoe',9999)
93      if(.not.ma_push_get(mt_dbl,5*mnoe,'dnoe',l_dnoe,i_dnoe))
94     + call md_abort('Memory allocation failed for dnoe',9999)
95      if(.not.ma_push_get(mt_int,mat,'anum',l_num,i_num))
96     + call md_abort('Memory allocation failed for anum',9999)
97      if(.not.ma_push_get(mt_dbl,mat,'amass',l_mas,i_mas))
98     + call md_abort('Memory allocation failed for amass',9999)
99      if(.not.ma_push_get(mt_byte,16*mwa,'cwa',l_cwa,i_cwa))
100     + call md_abort('Memory allocation failed for cwa',9999)
101      if(.not.ma_push_get(mt_dbl,mwa,'qwa',l_qwa,i_qwa))
102     + call md_abort('Memory allocation failed for qwa',9999)
103      if(.not.ma_push_get(mt_byte,16*msa,'csa',l_csa,i_csa))
104     + call md_abort('Memory allocation failed for csa',9999)
105      if(.not.ma_push_get(mt_int,mwa,'wat',l_wat,i_wat))
106     + call md_abort('Memory allocation failed for wat',9999)
107      if(.not.ma_push_get(mt_int,msa,'sat',l_sat,i_sat))
108     + call md_abort('Memory allocation failed for sat',9999)
109      if(.not.ma_push_get(mt_int,msa,'sgm',l_sgm,i_sgm))
110     + call md_abort('Memory allocation failed for sgm',9999)
111      if(.not.ma_push_get(mt_int,msa,'sfr',l_sfr,i_sfr))
112     + call md_abort('Memory allocation failed for sfr',9999)
113      if(.not.ma_push_get(mt_int,msa,'sml',l_sml,i_sml))
114     + call md_abort('Memory allocation failed for sml',9999)
115      if(.not.ma_push_get(mt_dbl,msa,'qsa',l_qsa,i_qsa))
116     + call md_abort('Memory allocation failed for qsa',9999)
117      if(.not.ma_push_get(mt_int,2*msb,'idsb',l_isb,i_isb))
118     + call md_abort('Memory allocation failed for idsb',9999)
119      if(.not.ma_push_get(mt_dbl,msb,'cdsb',l_csb,i_csb))
120     + call md_abort('Memory allocation failed for cdsb',9999)
121      if(.not.ma_push_get(mt_int,mseq,'lseq',l_lseq,i_lseq))
122     + call md_abort('Memory allocation failed for lseq',9999)
123      if(.not.ma_push_get(mt_int,msa,'ihop',l_ihop,i_ihop))
124     + call md_abort('Memory allocation failed for ihop',9999)
125      if(.not.ma_push_get(mt_int,msa,'istat',l_istat,i_istat))
126     + call md_abort('Memory allocation failed for istat',9999)
127c
128      nwa=0
129      nsa=0
130      nsm=0
131      noe=0
132c
133      if(util_print('where',print_debug)) then
134      write(lfnout,1002)
135 1002 format(' Memory allocated')
136      endif
137c
138c     read topology file
139c     ------------------
140c
141      if(.not.argos_prep_toprd(lfntop,filtop,lfnout,
142     + int_mb(i_num),dbl_mb(i_mas),mat,
143     + byte_mb(i_cwa),dbl_mb(i_qwa),mwa,nwa,
144     + int_mb(i_wat),int_mb(i_sat),int_mb(i_sgm),int_mb(i_sml),
145     + int_mb(i_sfr),
146     + byte_mb(i_csa),dbl_mb(i_qsa),msa,nsa,nsm,int_mb(i_isb),
147     + dbl_mb(i_csb),msb,nsb,dbl_mb(i_qu),nqu,slvnam,
148     + mseq,nseq,int_mb(i_lseq),int_mb(i_ihop),int_mb(i_istat)))
149     + call md_abort('argos_prep_toprd failed',9999)
150      if(util_print('coordinates',print_default)) then
151      write(lfnout,1003) filtop(1:index(filtop,' ')-1)
152 1003 format(' Topology',t40,a,/)
153      endif
154c
155      msm=nsm
156c
157      if(.not.ma_push_get(mt_int,4*msm,'eq',l_eq,i_eq))
158     + call md_abort('Memory allocation failed for eq',9999)
159c
160c     determine equivalent solute molecules
161c     -------------------------------------
162c
163      if(.not.argos_prep_equiv(byte_mb(i_csa),int_mb(i_sml),
164     + int_mb(i_sgm),
165     + msa,nsa,int_mb(i_eq),msm,nsm))
166     + call md_abort('Solute molecule equivalent test failed',0)
167c
168c     open the command file
169c     ---------------------
170c
171      open(unit=lfncmd,file=filcmd(1:index(filcmd,' ')-1),
172     + form='formatted',status='old',err=999)
173      rewind(lfncmd)
174c
175    1 continue
176c
177      read(lfncmd,1101,end=9,err=9997) card
178 1101 format(a)
179c
180      if(util_print('where',print_high)) then
181      write(lfnout,1099) card
182 1099 format('command card ',a)
183      endif
184c
185c     read distance restraints
186c     ========================
187c
188      if(card(1:6).eq.'disres') then
189      read(lfncmd,1101,end=9,err=9997) card2
190      if(.not.argos_prep_disres(card,card2,int_mb(i_sgm),
191     + byte_mb(i_csa),msa,
192     + nsa,int_mb(i_inoe),dbl_mb(i_dnoe),mnoe,nnoe))
193     + call md_abort('argos_prep_disres failed',9999)
194      endif
195c
196c     read pmf replication
197c     ====================
198c
199      if(card(1:7).eq.'rep_pmf') then
200      read(card(8:22),'(i5,i10)') npmfr,nsarep
201      ipmfr=numpmf
202      endif
203c
204c     read group selection
205c     ====================
206c
207      if(card(1:6).eq.'select') then
208      if(.not.argos_prep_select(card,int_mb(i_sgm),int_mb(i_sml),
209     + byte_mb(i_csa),msa,
210     + nsa,maxgrp,maxatm,igroup,lgroup,mgroup))
211     + call md_abort('argos_prep_select failed',9999)
212      endif
213c
214c     read pmf
215c     ========
216c
217      if(card(1:3).eq.'pmf') then
218      numpmf=numpmf+1
219      if(numpmf.gt.maxpmf) call md_abort('Increase maxpmf',maxpmf)
220      ipmf(numpmf,1)=0
221      if(card(5:7).eq.'con') ipmf(numpmf,1)=-1
222      if(card(5:7).eq.'dis') ipmf(numpmf,1)=1
223      if(card(5:7).eq.'ang') ipmf(numpmf,1)=2
224      if(card(5:7).eq.'tor') ipmf(numpmf,1)=3
225      if(card(5:7).eq.'imp') ipmf(numpmf,1)=4
226      if(card(5:7).eq.'lin') ipmf(numpmf,1)=5
227      if(card(5:7).eq.'pla') ipmf(numpmf,1)=6
228      if(card(5:7).eq.'bas') ipmf(numpmf,1)=7
229      if(card(5:7).eq.'zax') ipmf(numpmf,1)=8
230      if(card(5:7).eq.'zdi') ipmf(numpmf,1)=9
231      if(card(5:7).eq.'ZAX') ipmf(numpmf,1)=10
232      read(card(8:90),'(7i5,2f12.6,2e12.5)')
233     + (iopt(numpmf,i),i=1,3),(ipmf(numpmf,i),i=2,5),
234     + (rpmf(numpmf,i),i=1,4)
235      endif
236c
237      if(card(1:3).ne.'end') goto 1
238c
239    9 continue
240      close(unit=lfncmd,err=999)
241c
242      if(nnoe.gt.0) then
243      call argos_prep_wrtnoe(lfntop,filtop,mnoe,nnoe,
244     + int_mb(i_inoe),dbl_mb(i_dnoe))
245      write(lfnout,1004)
246 1004 format(/,' Distance restraints appended to topology')
247      endif
248c
249      if(numpmf.gt.0) then
250      if(.not.ma_push_get(mt_int,msa,'temp',l_temp,i_temp))
251     + call md_abort('Memory allocation failed for temp',9999)
252      if(.not.ma_push_get(mt_int,msa,'tempj',l_tempj,i_tempj))
253     + call md_abort('Memory allocation failed for tempj',9999)
254      call argos_prep_wrtpmf(lfntop,filtop,numpmf,
255     + maxpmf,ipmf,rpmf,iopt,maxgrp,maxatm,igroup,lgroup,mgroup,npmfr,
256     + ipmfr,nsarep,msa,nsa,int_mb(i_sgm),int_mb(i_sml),byte_mb(i_csa),
257     + int_mb(i_eq),msm,nsm,int_mb(i_temp),int_mb(i_tempj))
258      if(.not.ma_pop_stack(l_tempj))
259     + call md_abort('Memory deallocation failed for tempj',9999)
260      if(.not.ma_pop_stack(l_temp))
261     + call md_abort('Memory deallocation failed for temp',9999)
262      write(lfnout,1005)
263 1005 format(/,' Potential of mean force appended to topology')
264      endif
265c
266  999 continue
267c
268c     deallocate memory
269c     -----------------
270c
271      if(.not.ma_pop_stack(l_eq))
272     + call md_abort('Memory deallocation failed for eq',9999)
273      if(.not.ma_pop_stack(l_istat))
274     + call md_abort('Memory deallocation failed for istat',9999)
275      if(.not.ma_pop_stack(l_ihop))
276     + call md_abort('Memory deallocation failed for ihop',9999)
277      if(.not.ma_pop_stack(l_lseq))
278     + call md_abort('Memory deallocation failed for lseq',9999)
279      if(.not.ma_pop_stack(l_csb))
280     + call md_abort('Memory deallocation failed for cdsb',9999)
281      if(.not.ma_pop_stack(l_isb))
282     + call md_abort('Memory deallocation failed for idsb',9999)
283      if(.not.ma_pop_stack(l_qsa))
284     + call md_abort('Memory deallocation failed for qsa',9999)
285      if(.not.ma_pop_stack(l_sml))
286     + call md_abort('Memory deallocation failed for qsa',9999)
287      if(.not.ma_pop_stack(l_sfr))
288     + call md_abort('Memory deallocation failed for sfr',9999)
289      if(.not.ma_pop_stack(l_sgm))
290     + call md_abort('Memory deallocation failed for sgm',9999)
291      if(.not.ma_pop_stack(l_sat))
292     + call md_abort('Memory deallocation failed for sat',9999)
293      if(.not.ma_pop_stack(l_wat))
294     + call md_abort('Memory deallocation failed for wat',9999)
295      if(.not.ma_pop_stack(l_csa))
296     + call md_abort('Memory deallocation failed for csa',9999)
297      if(.not.ma_pop_stack(l_qwa))
298     + call md_abort('Memory deallocation failed for qwa',9999)
299      if(.not.ma_pop_stack(l_cwa))
300     + call md_abort('Memory deallocation failed for cwa',9999)
301      if(.not.ma_pop_stack(l_mas))
302     + call md_abort('Memory deallocation failed for amass',9999)
303      if(.not.ma_pop_stack(l_num))
304     + call md_abort('Memory deallocation failed for anum',9999)
305      if(.not.ma_pop_stack(l_dnoe))
306     + call md_abort('Memory deallocation failed for dnoe',9999)
307      if(.not.ma_pop_stack(l_inoe))
308     + call md_abort('Memory deallocation failed for inoe',9999)
309      if(.not.ma_pop_stack(l_qu))
310     + call md_abort('Memory deallocation failed for qu',9999)
311c
312      argos_prep_mknoe=.true.
313      return
314c
315 9997 continue
316      write(lfnout,1017) filcmd(1:index(filcmd,' ')-1)
317 1017 format(' Error reading commands',t40,a,/)
318      argos_prep_mknoe=.false.
319      return
320c
321 9999 continue
322      argos_prep_mknoe=.false.
323      return
324      end
325