1      subroutine qmmm_bq_data_init(irtdb)
2      implicit none
3c
4#include "mafdecls.fh"
5#include "rtdb.fh"
6#include "bq.fh"
7#include "errquit.fh"
8#include "qmmm_bq_data.fh"
9      integer irtdb
10c
11      character*30 bq_expand
12      character*30 pname
13
14      pname = "qmmm_bq_data_init"
15      nbqs=0
16      nbqw=0
17      nbq=nbqs+nbqw
18
19       if (.not.rtdb_get(irtdb,'qmmm:bq_exclude',mt_int,1,bq_exclude))
20     + call errquit(pname,bq_exclude,
21     &       RTDB_ERR)
22
23      if (.not. rtdb_get(irtdb,'qmmm:bq_update',mt_int,1,bq_update))
24     +    bq_update = 0
25
26      readbq = .false.
27      if (rtdb_cget(irtdb,"qmmm:bq:load:file",1,bqfile_in)) then
28        if(bqfile_in.ne."none") then
29         readbq = .true.
30         if(bqfile_in.eq."default") then
31           call util_file_name("bq.ind",.false.,.false.,bqfile_in)
32         else
33           call util_file_name_resolve(bqfile_in,.false.)
34         end if
35        end if
36      end if
37
38
39      bqsave = .false.
40      if (rtdb_cget(irtdb,'qmmm:bq:save:file',1,bqfile_out)) then
41        if(bqfile_out.ne."none") then
42         bqsave = .true.
43         if(bqfile_out.eq."default") then
44           call util_file_name("bq.ind",.false.,.false.,bqfile_out)
45         else
46           call util_file_name_resolve(bqfile_out,.false.)
47         end if
48        end if
49      end if
50c
51      if (.not.rtdb_cget(irtdb,"qmmm:bq_expand",1,bq_expand))
52     >      bq_expand = "none"
53      if(bq_expand.eq."none") then
54        allbqs = .false.
55        allbqw = .false.
56      else if(bq_expand.eq."all") then
57        allbqs = .true.
58        allbqw = .true.
59      else if(bq_expand.eq."solute") then
60        allbqs = .true.
61        allbqw = .false.
62      else if(bq_expand.eq."solvent") then
63        allbqs = .false.
64        allbqw = .true.
65      else
66       call errquit(pname,"unknown option for charge expand",0,0)
67      end if
68c
69c      if (.not.rtdb_get(irtdb,'qmmm:bq-solute',mt_log,1,allbqs))
70c     >     allbqs = .false.
71cc
72c      if (.not.rtdb_get(irtdb,'qmmm:bq-solvent',mt_log,1,allbqw))
73c     >     allbqw = .false.
74
75      if(.not.bq_create("qmmm charges",bq_handle))
76     + call errquit(pname//'Failed bq_create',0,CALC_ERR)
77
78      call qmmm_bq_data_load()
79c
80c    force separate calculation of bq energy
81c    ---------------------------------------
82      if (.not. rtdb_put(irtdb,'dft:bq_energy',mt_dbl,1,0.0d0))
83     + call errquit(pname//'setting dft:bq_energy',0,CALC_ERR)
84c
85
86      end
87
88      subroutine qmmm_bq_data_load()
89      implicit none
90c
91#include "mafdecls.fh"
92#include "errquit.fh"
93#include "qmmm_bq_data.fh"
94#include "qmmm.fh"
95#include "qmmm_utils.fh"
96#include "qmmm_params.fh"
97#include "global.fh"
98#include "bq.fh"
99#include "util.fh"
100#include "rtdb.fh"
101      integer i,k
102      integer i_lb
103      integer i_mbq,h_mbq
104      integer nlink
105      integer il,j
106      double precision charge
107      integer mwa
108      character*32 pname
109
110      pname = "qmmm_bq_data_load"
111
112      if(qmmm_print_debug())
113     >       write(*,*) "in",pname
114
115      call qmmm_bq_data_dealloc()
116c
117c     reset bq flags for solute and solvent
118c     -------------------------------------
119      if(allbqs) call mm_activate_bqszone()
120      if(allbqw) call mm_activate_bqwzone()
121c
122c     exclude bq's as specified by qmmm input
123c     ---------------------------------------
124      call mm_prune_bqzone()
125c
126c     get numbers of solute and solvent bq's
127c     --------------------------------------
128      if(readbq) then
129         call qmmm_read_nbq()
130      else
131         call mm_get_tot_nbqs(nbqs)
132         call mm_get_tot_nbqw(nbqw)
133      end if
134      call mm_get_mwa(mwa)
135      nbq=nbqs+nbqw
136      mbqw = nbqw/mwa
137
138      if(nbq.eq.0) return
139
140c     allocate bq structure
141c     ---------------------
142      call qmmm_bq_data_alloc()
143
144c     get global index for bq charges
145c     -------------------------------
146      if(readbq) then
147        call qmmm_read_ibq()
148        call mm_reset_solute_bqzone(nbqs,int_mb(i_ibq))
149        call mm_reset_solvent_bqzone(nbqw,int_mb(i_ibq+nbqs))
150      else
151        if(nbqs.ne.0) then
152          call mm_get_solute_bq_ind(nbqs,int_mb(i_ibq))
153        end if
154        if(nbqw.ne.0) then
155          call mm_get_solvent_ind_bq(nbqw,int_mb(i_ibq+nbqs))
156        end if
157      end if
158c
159c     get solute bq coord and charges from MM
160c     ---------------------------------------
161      if(nbqs.ne.0) then
162        call mm_get_solute_geom_bq(nbqs,
163     >                       int_mb(i_ibq),
164     >                       dbl_mb(i_cbq),
165     >                       dbl_mb(i_qbq))
166
167      end if
168
169      if(qmmm_print_debug())
170     >       write(*,*) "finished setting up solute bq charges"
171c
172c     append solvent bq coord and charges from MM
173c     -------------------------------------------
174      if(nbqw.ne.0) then
175        call mm_get_solvent_geom_bq(nbqw,int_mb(i_ibq+nbqs),
176     >                       dbl_mb(i_cbq+3*nbqs),
177     >                       dbl_mb(i_qbq+nbqs))
178      end if
179
180      if(qmmm_print_debug())
181     >       write(*,*) "finished setting up solvent bq charges"
182
183c
184      if(qmmm_master()) then
185         write(*,*) "Total number of Bq charges",nbq
186         write(*,*) "Number of solute Bq charges",nbqs
187         write(*,*) "Number of solvent Bq charges",nbqw
188      end if
189c
190c     save bq index into the file
191c     ---------------------------
192      if(bqsave)
193     +  call qmmm_save_bq_index()
194c
195c     create and activate bq structure in bq module
196c     ---------------------------------------------
197      if(.not.bq_pset(bq_handle,nbq,h_qbq,h_cbq))
198     + call errquit(pname//'Failed bq_pset',0,CALC_ERR)
199      if(.not.bq_activate(bq_handle))
200     + call errquit(pname//'Failed bq_activate',0,CALC_ERR)
201      charge = 0.0d0
202      do i=1,nbq
203       charge = charge + dbl_mb(i_qbq+i-1)
204      end do
205      if(ga_nodeid().eq.0) then
206       write(*,*) "Total Bq charge: ", charge
207      end if
208
209      if(util_print('bq_geom', print_high))
210     + call bq_print_info(bq_handle)
211
212       if(ga_nodeid().eq.0)
213     > call qmmm_print_pdb_bq(nbq,"-bq.pdb",
214     >                    dbl_mb(i_cbq),
215     >                    dbl_mb(i_qbq))
216
217      if(qmmm_print_debug()) then
218        call util_print_centered(6,"Bq information",
219     >              36, .true.)
220        write(*,'(/,A,T20,A,/)') "global index","coordinates(au)"
221        do i=1,nbqs+nbqw
222         write(*,'(I5,T20,3F12.6)') int_mb(i_ibq+i-1),
223     >              (dbl_mb(i_cbq+3*(i-1)+k-1),k=1,3)
224        end do
225      end if
226
227      if(qmmm_print_debug())
228     >       write(*,*) "out",pname
229
230      end
231
232      subroutine qmmm_bq_data_reload()
233      implicit none
234c
235#include "mafdecls.fh"
236#include "errquit.fh"
237#include "qmmm_bq_data.fh"
238#include "qmmm.fh"
239      integer counter
240      character*32 pname
241      save    counter
242      DATA counter /0/
243
244      pname = "qmmm_bq_data_reload"
245
246      if(qmmm_print_debug())
247     >       write(*,*) "in",pname
248
249      counter = counter + 1
250      if(counter.eq.bq_update) then
251        call  qmmm_bq_data_load()
252        counter = 0
253      end if
254
255      if(qmmm_print_debug())
256     >       write(*,*) "out",pname
257
258      return
259      end
260
261      subroutine qmmm_bq_data_update_active()
262      implicit none
263c
264#include "mafdecls.fh"
265#include "errquit.fh"
266#include "qmmm_bq_data.fh"
267#include "qmmm.fh"
268#include "qmmm_utils.fh"
269#include "qmmm_params.fh"
270#include "global.fh"
271#include "bq.fh"
272#include "util.fh"
273      integer i,k
274      character*32 pname
275
276      pname = "qmmm_bq_data_create_active"
277
278      if(qmmm_print_debug())
279     >       write(*,*) "in",pname
280c
281      do i=1,nbq
282         log_mb(i_abq+i-1)= .false.
283      end do
284c
285c     get active list of bqs
286c     ------------------------
287      if(nbqs.ne.0)
288     >   call qmmm_cons_get_actmaps(nbqs,int_mb(i_ibq),
289     >                       log_mb(i_abq))
290
291      if(nbqw.ne.0)
292     >  call  qmmm_cons_get_actmapw(nbqw,int_mb(i_ibq+nbqs),
293     >                       log_mb(i_abq+nbqs))
294
295      nbqa = 0
296      do i=1,nbq
297         if(log_mb(i_abq+i-1)) then
298           nbqa = nbqa + 1
299         end if
300      end do
301
302      if(ga_nodeid().eq.0) then
303       write(*,*) "Total number of active Bq charges ",nbqa
304      end if
305
306      if(qmmm_print_debug())
307     >       write(*,*) "out",pname
308
309      end
310
311      subroutine qmmm_bq_data_reset_active()
312      implicit none
313c
314#include "mafdecls.fh"
315#include "errquit.fh"
316#include "qmmm_bq_data.fh"
317#include "qmmm.fh"
318#include "qmmm_utils.fh"
319#include "qmmm_params.fh"
320#include "global.fh"
321#include "bq.fh"
322#include "util.fh"
323      integer i,k
324      character*32 pname
325
326      pname = "qmmm_bq_data_create_active"
327
328      if(qmmm_print_debug())
329     >       write(*,*) "in",pname
330c
331      do i=1,nbq
332         log_mb(i_abq+i-1)= .true.
333      end do
334      nbqa = nbq
335      do i=1,nbq
336         if(log_mb(i_abq+i-1)) then
337           nbqa = nbqa + 1
338         end if
339      end do
340
341      if(ga_nodeid().eq.0) then
342       write(*,*) "number of active bqs ",nbqa
343      end if
344
345      if(qmmm_print_debug())
346     >       write(*,*) "out",pname
347
348      end
349
350      subroutine qmmm_bq_print_info()
351      implicit none
352c
353#include "mafdecls.fh"
354#include "errquit.fh"
355#include "qmmm_bq_data.fh"
356#include "qmmm.fh"
357#include "qmmm_utils.fh"
358#include "qmmm_params.fh"
359#include "global.fh"
360#include "bq.fh"
361#include "util.fh"
362      integer i,k
363      character*32 pname
364
365      pname = "qmmm_bq_print_info"
366
367      if(qmmm_print_debug())
368     >       write(*,*) "in",pname
369
370        call util_print_centered(6,"Bq information",
371     >              36, .true.)
372        write(*,'(/,A,T20,A,/)') "global index","coordinates(au)"
373        do i=1,nbqs+nbqw
374         write(*,'(I5,T20,3F12.6)') int_mb(i_ibq+i-1),
375     >              (dbl_mb(i_cbq+3*(i-1)+k-1),k=1,3)
376        end do
377
378      if(qmmm_print_debug())
379     >       write(*,*) "out",pname
380
381      end
382
383      subroutine qmmm_bq_coord_update()
384      implicit none
385c
386#include "mafdecls.fh"
387#include "errquit.fh"
388#include "qmmm_bq_data.fh"
389#include "qmmm.fh"
390#include "global.fh"
391c
392      integer i
393      if(nbqs.ne.0) then
394        call mm_get_solute_coord(nbqs,
395     >                       int_mb(i_ibq),
396     >                       dbl_mb(i_cbq))
397
398      end if
399
400c
401c     append solvent bq coord and charges from MM
402c     -------------------------------------------
403      if(nbqw.ne.0) then
404        call mm_get_solvent_coord(nbqw,int_mb(i_ibq+nbqs),
405     >                       dbl_mb(i_cbq+3*nbqs))
406      end if
407
408      if(qmmm_print_debug())
409     >       write(*,*) "finshed updating up bq coords"
410
411      end
412
413      subroutine qmmm_bq_data_alloc()
414      implicit none
415c
416#include "mafdecls.fh"
417#include "errquit.fh"
418#include "qmmm_bq_data.fh"
419      integer i
420      if(nbq.eq.0) return
421
422      if(.not.ma_alloc_get(mt_log,nbq,'abq',h_abq,i_abq))
423     + call errquit('qmmm: Failed to allocate memory for abq',nbq,
424     +       MA_ERR)
425      do i=1,nbq
426        log_mb(i_abq+i-1)=.true.
427      end do
428
429      if(.not.ma_alloc_get(mt_dbl,3*nbq,'cbq',h_cbq,i_cbq))
430     + call errquit('qmmm: Failed to allocate memory for cbq',3*nbq,
431     +       MA_ERR)
432      call dfill(3*nbq,0.d0,dbl_mb(i_cbq),1)
433      if(.not.ma_alloc_get(mt_dbl,nbq,'qbq',h_qbq,i_qbq))
434     + call errquit('qmmm: Failed to allocate memory for qbq',nbq,
435     &       MA_ERR)
436      call dfill(nbq,0.d0,dbl_mb(i_qbq),1)
437      if(.not.ma_alloc_get(mt_int,nbq,'ibq',h_ibq,i_ibq))
438     + call errquit('qmmm: Failed to allocate memory for ibq',nbq,
439     &       MA_ERR)
440      call ifill(nbq,0,int_mb(i_ibq),1)
441
442      if(.not.ma_alloc_get(mt_int,nbq,'i_ibqa',h_ibqa,i_ibqa))
443     + call errquit('qmmm: Failed to allocate memory for i_ibqa',nbq,
444     &       MA_ERR)
445      call ifill(nbq,0,int_mb(i_ibqa),1)
446
447      end
448
449      subroutine qmmm_bq_data_dealloc()
450      implicit none
451c
452#include "mafdecls.fh"
453#include "errquit.fh"
454#include "qmmm_bq_data.fh"
455#include "bq.fh"
456      logical ignore
457
458      if(nbq.eq.0) return
459
460      if(.not.ma_free_heap(h_ibqa))
461     + call errquit('qmmm:memory deallocation for ibqa',nbqa,
462     &       MA_ERR)
463
464      if(.not.ma_free_heap(h_ibq))
465     + call errquit('qmmm:memory deallocation for ibq',nbq,
466     &       MA_ERR)
467
468      if(.not.ma_free_heap(h_qbq))
469     + call errquit('qmmm:memory deallocation for qbq',nbq,
470     &       MA_ERR)
471
472      if(.not.ma_free_heap(h_cbq))
473     + call errquit('qmmm: memory deallocation for cbq',3*nbq,
474     +       MA_ERR)
475
476      if(.not.ma_free_heap(h_abq))
477     + call errquit('qmmm: memory deallocation for cbq',3*nbq,
478     +       MA_ERR)
479
480       nbqs=0
481       nbqw=0
482       nbq=0
483
484       ignore = bq_destroy(bq_handle)
485
486
487      end
488
489      function qmmm_get_nbq()
490      implicit none
491#include "qmmm_bq_data.fh"
492
493      integer qmmm_get_nbq
494
495      qmmm_get_nbq = nbqs+nbqw
496
497      end
498
499      function qmmm_get_nbqa()
500      implicit none
501#include "qmmm_bq_data.fh"
502
503      integer qmmm_get_nbqa
504
505      qmmm_get_nbqa = nbqa
506
507      end
508
509      function qmmm_get_nbqs()
510      implicit none
511#include "qmmm_bq_data.fh"
512
513      integer qmmm_get_nbqs
514
515      qmmm_get_nbqs = nbqs
516
517      end
518
519      function qmmm_get_nbqw()
520      implicit none
521#include "qmmm_bq_data.fh"
522
523      integer qmmm_get_nbqw
524
525      qmmm_get_nbqw = nbqw
526
527      end
528
529      function qmmm_get_h_qbq()
530      implicit none
531#include "qmmm_bq_data.fh"
532      integer qmmm_get_h_qbq
533
534      qmmm_get_h_qbq = h_qbq
535
536      end
537
538      function qmmm_get_i_qbq()
539      implicit none
540#include "qmmm_bq_data.fh"
541      integer qmmm_get_i_qbq
542
543      qmmm_get_i_qbq = i_qbq
544
545      end
546
547      function qmmm_get_h_cbq()
548      implicit none
549#include "qmmm_bq_data.fh"
550      integer qmmm_get_h_cbq
551
552      qmmm_get_h_cbq = h_cbq
553
554      end
555
556      function qmmm_get_i_cbq()
557      implicit none
558#include "qmmm_bq_data.fh"
559      integer qmmm_get_i_cbq
560
561      qmmm_get_i_cbq = i_cbq
562
563      end
564
565      function qmmm_get_i_ibq()
566      implicit none
567#include "qmmm_bq_data.fh"
568      integer qmmm_get_i_ibq
569
570      qmmm_get_i_ibq = i_ibq
571
572      end
573
574      function qmmm_get_i_ibqa()
575      implicit none
576#include "qmmm_bq_data.fh"
577      integer qmmm_get_i_ibqa
578
579      qmmm_get_i_ibqa = i_ibqa
580
581      end
582
583      function qmmm_get_h_cbqs()
584      implicit none
585#include "qmmm_bq_data.fh"
586      integer qmmm_get_h_cbqs
587
588      qmmm_get_h_cbqs = h_cbqs
589
590      end
591
592      function qmmm_get_i_cbqs()
593      implicit none
594#include "qmmm_bq_data.fh"
595      integer qmmm_get_i_cbqs
596
597      qmmm_get_i_cbqs = i_cbqs
598
599      end
600
601      function qmmm_get_h_cbqw()
602      implicit none
603#include "qmmm_bq_data.fh"
604      integer qmmm_get_h_cbqw
605
606      qmmm_get_h_cbqw = h_cbqw
607
608      end
609
610      function qmmm_get_i_cbqw()
611      implicit none
612#include "qmmm_bq_data.fh"
613      integer qmmm_get_i_cbqw
614
615      qmmm_get_i_cbqw = i_cbqw
616
617      end
618
619      function qmmm_get_i_qbqw()
620      implicit none
621#include "qmmm_bq_data.fh"
622      integer qmmm_get_i_qbqw
623
624      qmmm_get_i_qbqw = i_qbqw
625
626      end
627
628      function qmmm_get_bq_exclude()
629      implicit none
630#include "qmmm_bq_data.fh"
631
632      integer qmmm_get_bq_exclude
633
634      qmmm_get_bq_exclude = bq_exclude
635
636      end
637
638      subroutine qmmm_read_ibq()
639      implicit none
640#include "stdio.fh"
641#include "errquit.fh"
642#include "qmmm.fh"
643#include "util.fh"
644#include "inp.fh"
645#include "mafdecls.fh"
646#include "global.fh"
647#include "msgids.fh"
648#include "qmmm_bq_data.fh"
649
650c     local variables
651      integer nf
652      character*30 pname
653
654      integer i
655      integer fn
656      integer anbq,anbqs,anbqw
657
658      pname = "qmmm_read_ibq"
659
660      if(qmmm_master()) then
661        if(.not.qmmm_get_io_unit(fn))
662     >       call errquit("cannot get file number",0,0)
663        open(unit=fn,status="old",form="formatted",file=bqfile_in)
664        read(fn,*) anbqs,anbqw,anbq
665        do i=1,anbq
666          read(fn,*) int_mb(i_ibq+i-1)
667        end do
668        close(fn)
669        call util_print_centered(luout,
670     >              "Bq index has been read from external file "//
671     >               bqfile_in,
672     >              36, .true.)
673      end if
674      call ga_brdcst(msg_qmmm_ind,int_mb(i_ibq),
675     >      nbq*ma_sizeof(mt_int,1,mt_byte),0)
676
677      call ga_sync()
678
679      end
680
681      subroutine qmmm_read_nbq()
682      implicit none
683#include "errquit.fh"
684#include "qmmm.fh"
685#include "util.fh"
686#include "inp.fh"
687#include "mafdecls.fh"
688#include "global.fh"
689#include "msgids.fh"
690#include "qmmm_bq_data.fh"
691
692c     local variables
693      integer nf
694
695      integer i
696      integer fn
697
698      if(qmmm_master()) then
699        if(.not.qmmm_get_io_unit(fn))
700     >       call errquit("cannot get file number",0,0)
701        open(unit=fn,status="old",form="formatted",file=bqfile_in)
702        read(fn,*) nbqs,nbqw
703        close(fn)
704      end if
705
706      call ga_brdcst(msg_qmmm_nbqs,nbqs,ma_sizeof(mt_int,1,mt_byte),0)
707      call ga_brdcst(msg_qmmm_nbqw,nbqw,ma_sizeof(mt_int,1,mt_byte),0)
708
709
710      end
711
712      subroutine qmmm_save_bq_index()
713      implicit none
714#include "stdio.fh"
715#include "errquit.fh"
716#include "qmmm.fh"
717#include "util.fh"
718#include "global.fh"
719#include "inp.fh"
720#include "qmmm_bq_data.fh"
721#include "mafdecls.fh"
722
723c     local variables
724      integer nf
725      character*(nw_max_path_len) filename
726
727      integer i
728      integer fn
729
730      if(qmmm_master()) then
731
732      if(.not.qmmm_get_io_unit(fn))
733     >       call errquit("cannot get file number",0,0)
734
735      open(unit=fn,status="unknown",form="formatted",file=bqfile_out)
736
737      write(fn,*) nbqs,nbqw,nbq
738      do i=1,nbq
739       write(fn,*) int_mb(i_ibq+i-1)
740      end do
741
742      call util_flush(fn)
743      close(fn)
744
745      call util_print_centered(luout,
746     >              "Bq index has been saved to external file "//
747     >               bqfile_out,
748     >              36, .true.)
749
750      end if
751
752      call ga_sync()
753
754      end
755
756      subroutine qmmm_bq_push_grad(nt,
757     >                       g)
758      implicit none
759
760#include "qmmm_bq_data.fh"
761#include "qmmm_params.fh"
762#include "global.fh"
763#include "mafdecls.fh"
764#include "errquit.fh"
765
766      integer nt
767      double precision g(3*nt)
768c
769      integer i_ibqg,h_ibqg
770      integer i,j,k,ia
771      integer ii,jj
772      integer nbqsa,nbqwa,mwa
773      character*32 pname
774
775      pname = "qmmm_bq_push_grad"
776      if(.not.ma_push_get(mt_int,nbqa,'ibqg',h_ibqg,i_ibqg))
777     +   call errquit( pname//'Failed to allocate memory for ibqs',
778     +   nbqa, MA_ERR)
779
780
781      nbqsa = 0
782      do i=1,nbqs
783         if(log_mb(i_abq+i-1)) then
784           nbqsa = nbqsa + 1
785           int_mb(i_ibqg+nbqsa-1) = int_mb(i_ibq+i-1)
786         end if
787      end do
788
789      call mm_get_mwa(mwa)
790      nbqwa = 0
791      do i=1,nbq-nbqs,mwa
792         if(log_mb(i_abq+nbqs+i-1)) then
793           nbqwa = nbqwa + 1
794           int_mb(i_ibqg+nbqsa+nbqwa-1) = int_mb(i_ibq+nbqs+(i-1)/mwa)
795         end if
796      end do
797
798
799      if(nbqsa.ne.0)
800     >  call mm_add_solute_force(nbqsa,
801     >                            int_mb(i_ibqg),
802     >                            g)
803
804
805      if(nbqwa.ne.0)
806     > call mm_add_solvent_force(mwa*nbqwa,
807     >                            int_mb(i_ibqg+nbqsa),
808     >                            g(3*nbqsa+1))
809
810
811      if(.not.ma_pop_stack(h_ibqg))
812     +   call errquit( pname//'Failed to deallocate memory for ibqg',
813     +   nbqa, MA_ERR)
814
815
816      end
817
818      subroutine qmmm_get_cqbqa(nt,
819     >                       c,
820     >                       q)
821      implicit none
822
823#include "qmmm_bq_data.fh"
824#include "global.fh"
825#include "mafdecls.fh"
826#include "errquit.fh"
827#include "qmmm_params.fh"
828
829      integer nt
830      double precision c(3,nt)
831      double precision q(nt)
832c
833      integer i,j,k,ia
834      integer ii,jj
835
836
837      ia = 0
838      do i=1,nbqs
839        if(log_mb(i_abq+i-1)) then
840           ia = ia + 1
841           q(ia) = dbl_mb(i_qbq+i-1)
842           do k=1,3
843             c(k,ia) = dbl_mb(i_cbq+(i-1)*3+k-1)
844           end do
845        end if
846      end do
847
848      do i=nbqs+1,nbq
849        if(log_mb(i_abq+i-1)) then
850           ia = ia + 1
851           q(ia) = dbl_mb(i_qbq+i-1)
852           do k=1,3
853             c(k,ia) = dbl_mb(i_cbq+(i-1)*3+k-1)
854           end do
855        end if
856      end do
857
858
859      end
860
861c $Id$
862