1
2      subroutine mm_geom_init(rtdb,geomname)
3      implicit none
4#include "mafdecls.fh"
5#include "errquit.fh"
6#include "geom.fh"
7#include "rtdb.fh"
8#include "util.fh"
9#include "global.fh"
10#include "inp.fh"
11#include "mm_coords_data.fh"
12#include "mm_bond_coords_data.fh"
13#include "mm_link_data.fh"
14#include "mm_geom_data.fh"
15
16      integer rtdb
17      character*(*) geomname
18
19      character*255 aname
20      character*32 pname
21      character*30 operation
22      logical orestrt
23      integer geom
24      integer i, j
25      logical ignore
26      integer i_c,h_c
27      integer i_tag,h_tag
28      integer i_m,h_m
29      integer i_atn,h_atn
30      integer i_dbq,h_dbq
31      double precision scale
32      logical geom_tag_to_atn
33      logical geom_tag_to_charge
34      external geom_tag_to_atn
35      external geom_tag_to_charge
36
37
38      pname = "mm_geom_init"
39c     write(*,*) pname
40
41      ignore = rtdb_cget(rtdb,'task:operation',1,operation)
42
43c     deallocate all previous allocated arrays just in case
44      call mm_geom_end()
45
46      nact = nqm
47      if(qmlink) nact = nqm + nlink
48
49      nfg = nqm + nlink
50
51      aux_geom = .false.
52c     TP: refer to qmmm.F for operation = "hessian" and
53c         operation = "freq". aux_geom = .true.
54      if(operation.eq."optimize")  aux_geom =.true.
55      if(operation.eq."neb")  aux_geom =.true.
56      if(operation.eq."freq")  aux_geom =.true.
57      if(operation.eq."hessian")  aux_geom =.true.
58      if(operation.eq."saddle")  aux_geom =.true.
59
60c     initialize indexing for a full qm geometry
61      if(.not.ma_alloc_get(mt_int,nact,'mm act atom ind',
62     &                     h_iact,i_iact))
63     & call errquit(pname//'Failed to allocate heap',nact,
64     &       MA_ERR)
65
66      if(.not.ma_alloc_get(mt_int,nfg,'mm fullg ind',h_ifg,i_ifg))
67     & call errquit(pname//'Failed to allocate heap',nfg,
68     &       MA_ERR)
69
70      call icopy(nqm,int_mb(i_iqm),1,int_mb(i_ifg),1)
71      call icopy(nlink,int_mb(i_lb+nlink),1,int_mb(i_ifg+nqm),1)
72
73      do i=1,nact
74        int_mb(i_iact+i-1) = int_mb(i_iqml+i-1)
75      end do
76
77      if(operation.eq."neb") then
78        if(.not.rtdb_get(rtdb,"mm:neb:restart",mt_log,1,orestrt))
79     >     orestrt = .false.
80        if(orestrt) return
81      end if
82
83
84      if(.not.ma_push_get(mt_dbl,3*nact,'c',h_c,i_c))
85     & call errquit('mm: Failed to allocate memory for c',
86     & 3*nact, MA_ERR)
87      if(.not.ma_push_get(mt_dbl,nact,'q',h_dbq,i_dbq))
88     & call errquit('mm: Failed to allocate memory for q',nact,
89     &       MA_ERR)
90      if(.not.ma_push_get(mt_dbl,nact,'m',h_m,i_m))
91     & call errquit('mm: Failed to allocate memory for m',nact,
92     &       MA_ERR)
93      if(.not.ma_push_get(mt_int,nact,'inum',h_atn,i_atn))
94     & call errquit('mm: Failed to allocate memory for atn',nact,
95     &       MA_ERR)
96      if(.not.ma_push_get(mt_byte,16*nact,'t',h_tag,i_tag))
97     & call errquit('mm: Failed to allocate memory for t',nact,
98     &       MA_ERR)
99
100c     assign coordinates for active atoms
101      do i=1,nact
102        dbl_mb(i_c+(i-1)*3)   = dbl_mb(i_rqml+3*(i-1))
103        dbl_mb(i_c+(i-1)*3+1) = dbl_mb(i_rqml+3*(i-1)+1)
104        dbl_mb(i_c+(i-1)*3+2) = dbl_mb(i_rqml+3*(i-1)+2)
105        call mm_coords_tag_set(byte_mb(i_tqml+16*(i-1)),
106     &                         byte_mb(i_tag+16*(i-1)))
107      end do
108
109      call util_convert_units("angstrom","au",scale)
110      call dscal(3*nact, scale,dbl_mb(i_c),1)
111
112      if(.not.geom_tag_to_charge(nact,byte_mb(i_tag),dbl_mb(i_dbq)))
113     & call errquit('mm: Failed to get charge from tag',nact,MA_ERR)
114
115      if(.not.geom_tag_to_atn(nact,byte_mb(i_tag),int_mb(i_atn)))
116     & call errquit('mm: Failed to get atn from tag',nact,MA_ERR)
117
118c     if nlink > 0 assign empirical charge to mmlink
119      if(qmlink) then
120        do i=1,nlink
121          dbl_mb(i_dbq+nqm+i-1) = dbl_mb(i_lnkchg+i-1)
122        end do
123      end if
124
125      call mm_atn_get_mass(nact,int_mb(i_atn),dbl_mb(i_m))
126
127c     ignore = rtdb_delete(rtdb,"geometry")
128      ignore = rtdb_delete(rtdb,geomname(1:inp_strlen(geomname)))
129
130c     if(.not.geom_create(geom,"geometry"))
131      if(.not.geom_create(geom,geomname(1:inp_strlen(geomname))))
132     & call errquit('mm: Failed to create geometry',0, GEOM_ERR)
133
134      if(.not.geom_cart_set(geom,nact,byte_mb(i_tag),
135     &                                dbl_mb(i_c),
136     &                                dbl_mb(i_dbq)))
137     & call errquit('mm: Failed to initialize geometry',0, GEOM_ERR)
138
139      if(.not.geom_masses_set(geom,nact,dbl_mb(i_m)))
140     & call errquit('mm: Failed to initialize masses',0, GEOM_ERR)
141      call geom_compute_values(geom)
142
143      if(.not.geom_rtdb_store(rtdb,geom,
144     &        geomname(1:inp_strlen(geomname))))
145     & call errquit('mm: Failed to store geom to rtdb',0, RTDB_ERR)
146
147      if(.not.geom_destroy(geom))
148     & call errquit('mm: Failed to destroy geometry',0, GEOM_ERR)
149
150      if(.not.ma_pop_stack(h_tag))
151     & call errquit('mm: Failed to deallocate stack t_all',nact,
152     &       MA_ERR)
153      if(.not.ma_pop_stack(h_atn))
154     & call errquit('mm: Failed to deallocate stack atn_all',nact,
155     &       MA_ERR)
156      if(.not.ma_pop_stack(h_m))
157     & call errquit('mm: Failed to deallocate stack m_all',nact,
158     &       MA_ERR)
159      if(.not.ma_pop_stack(h_dbq))
160     & call errquit('mm: Failed to deallocate stack q_all',nact,
161     &       MA_ERR)
162      if(.not.ma_pop_stack(h_c))
163     & call errquit('mm: Failed to deallocate stack c_all',nact,
164     &       MA_ERR)
165
166      end
167
168      subroutine mm_geom_end()
169      implicit none
170#include "mafdecls.fh"
171#include "errquit.fh"
172#include "geom.fh"
173#include "rtdb.fh"
174#include "global.fh"
175#include "inp.fh"
176#include "mm_geom_data.fh"
177
178      character*30 pname
179      pname = "mm_geom_end"
180
181      if(nfg.gt.0) then
182        if (.not.ma_free_heap(h_ifg))    goto 911
183        if (.not.ma_free_heap(h_iact))   goto 911
184        nfg = 0
185        nact = 0
186        aux_geom = .false.
187      end if
188
189      return
190911   call errquit("error "//trim(pname),0,-1)
191
192      end
193
194      subroutine mm_geom_create_full(rtdb)
195      implicit none
196#include "mafdecls.fh"
197#include "errquit.fh"
198#include "geom.fh"
199#include "rtdb.fh"
200#include "global.fh"
201#include "inp.fh"
202#include "mm_geom_data.fh"
203
204      character*32 pname
205      integer rtdb
206      integer i_c,h_c
207      integer i_t,h_t
208      integer i_m,h_m
209      integer i_atn,h_atn
210      integer i_dbq,h_dbq
211
212      pname = "geom_create_full"
213c     write(*,*) 'in ', pname
214
215      if(.not.ma_push_get(mt_dbl,3*nfg,'c',h_c,i_c))
216     & call errquit('mm: Failed to allocate memory for c',
217     & 3*nfg, MA_ERR)
218      if(.not.ma_push_get(mt_dbl,nfg,'q',h_dbq,i_dbq))
219     & call errquit('mm: Failed to allocate memory for q',nfg,
220     &       MA_ERR)
221      if(.not.ma_push_get(mt_dbl,nfg,'m',h_m,i_m))
222     & call errquit('mm: Failed to allocate memory for m',nfg,
223     &       MA_ERR)
224      if(.not.ma_push_get(mt_int,nfg,'inum',h_atn,i_atn))
225     & call errquit('mm: Failed to allocate memory for atn',nfg,
226     &       MA_ERR)
227      if(.not.ma_push_get(mt_byte,16*nfg,'t',h_t,i_t))
228     & call errquit('mm: Failed to allocate memory for t',nfg,
229     &       MA_ERR)
230
231      call mm_get_geom(rtdb,nfg,int_mb(i_ifg),dbl_mb(i_c),
232     &                 dbl_mb(i_dbq),dbl_mb(i_m),
233     &                 int_mb(i_atn),byte_mb(i_t))
234
235      if(.not.rtdb_cget(rtdb,'geometry',1,oldgeom))
236     & oldgeom = ' '
237
238      if(.not.rtdb_cput(rtdb,'geometry',1,'full_geom'))
239     & call errquit(pname//' storing geom name to rtdb',0, RTDB_ERR)
240
241c     release temporary memory
242c     ------------------------
243      if(.not.ma_pop_stack(h_t))
244     & call errquit('mm: Failed to deallocate stack t_all',nfg,
245     &       MA_ERR)
246      if(.not.ma_pop_stack(h_atn))
247     & call errquit('mm: Failed to deallocate stack atn_all',nfg,
248     &       MA_ERR)
249      if(.not.ma_pop_stack(h_m))
250     & call errquit('mm: Failed to deallocate stack m_all',nfg,
251     &       MA_ERR)
252      if(.not.ma_pop_stack(h_dbq))
253     & call errquit('mm: Failed to deallocate stack q_all',nfg,
254     &       MA_ERR)
255      if(.not.ma_pop_stack(h_c))
256     & call errquit('mm: Failed to deallocate stack c_all',nfg,
257     &       MA_ERR)
258
259
260      end
261
262      subroutine mm_get_geom(rtdb,nt,ifg,c,dbq,m,atn,t)
263      implicit none
264#include "mafdecls.fh"
265#include "errquit.fh"
266#include "geom.fh"
267#include "rtdb.fh"
268#include "util.fh"
269#include "mm_bond_coords_data.fh"
270#include "mm_coords_data.fh"
271
272      integer rtdb
273      integer geom
274      integer nt
275      integer ifg(nt)
276      double precision c(3,nt)
277      double precision dbq(nt)
278      double precision m(nt)
279      integer atn(nt)
280      character*16 t(nt)
281
282      integer i
283      integer nat
284      double precision scale
285      logical status
286      logical ignore
287      character*30 message
288      character*30 pname
289
290      pname = "mm_get_geom"
291
292c     if (.not. geom_create(geom, 'full_geom'))
293      if (.not. geom_create(geom, 'geometry'))
294     &     call errquit('cons_create_geom',0, GEOM_ERR)
295
296c     if (.not. geom_rtdb_load(rtdb, geom, 'full_geom'))
297      if (.not. geom_rtdb_load(rtdb, geom, 'geometry'))
298     &     call errquit('cons_load_geom',0, RTDB_ERR)
299
300      status=geom_ncent(geom,nat)
301
302      if(.not.status)
303     & call errquit('cons_init: geom_create?',70, GEOM_ERR)
304
305      if(.not.geom_cart_get2(geom,nat,t,c,dbq,atn))
306     &    goto 911
307
308      call util_convert_units("au","angstrom",scale)
309      call dscal(3*nat, scale,c,1)
310
311      if(.not.geom_destroy(geom))
312     &    goto 911
313
314      ignore = geom_rtdb_delete(rtdb,"full_geom")
315
316      if(.not.geom_create(geom,"full_geom"))
317     & call errquit('mm: Failed to create geometry',0, GEOM_ERR)
318
319      call mm_links_adjust(nt,ifg,atn,t,c,dbq)
320
321
322      call util_convert_units("angstrom","au",scale)
323      call dscal(3*nt, scale,c,1)
324
325      if(.not.geom_cart_set(geom,nt,t,c,dbq))
326     & call errquit('mm: Failed to initialize geometry',0, GEOM_ERR)
327
328      call mm_atn_get_mass(nt,atn,m)
329
330      if(.not.geom_masses_set(geom,nt,m))
331     & call errquit('mm: Failed to initialize masses',0, GEOM_ERR)
332      call geom_compute_values(geom)
333
334      if(.not.geom_rtdb_store(rtdb,geom,"full_geom"))
335     & call errquit('mm: Failed to store geom to rtdb',0, RTDB_ERR)
336
337      if(.not.geom_destroy(geom))
338     & call errquit('mm: Failed to destroy geometry',0, GEOM_ERR)
339
340      return
341911   call errquit("error "//trim(message),0,-1)
342      end
343
344      subroutine mm_atn_get_mass(n,atn,m)
345      implicit none
346#include "mafdecls.fh"
347#include "errquit.fh"
348#include "geom.fh"
349      integer n
350      integer atn(n)
351      double precision m(n)
352
353c     local variables
354      integer i
355      character*32 pname
356
357      pname = "mm_atn_get_mass"
358
359      do i=1,n
360        if(.not.geom_atn_to_default_mass(atn(i),m(i)))
361     &    call errquit(pname,0, GEOM_ERR)
362
363      end do
364
365      end
366
367      subroutine mm_geom_push_active(rtdb)
368      implicit none
369#include "mafdecls.fh"
370#include "errquit.fh"
371#include "geom.fh"
372#include "rtdb.fh"
373#include "global.fh"
374#include "inp.fh"
375#include "mm_geom_data.fh"
376#include "mm_bond_coords_data.fh"
377#include "mm_coords_data.fh"
378
379      integer rtdb
380
381      integer geom
382      integer ncent
383      integer i_ctmp,h_ctmp
384      integer i, j
385      integer nt
386      double precision scale
387      character*32 pname
388
389      pname = "mm_geom_push_active"
390c     write(*,*) 'in ', pname
391
392      nt = nact
393
394c     --------------------------------------
395c     get active coordinates out of geometry
396c     --------------------------------------
397      if(.not.geom_create(geom,'geometry'))
398     & call errquit('mm: Failed to create geometry',0, GEOM_ERR)
399
400      if(.not.geom_rtdb_load(rtdb,geom,"geometry"))
401     & call errquit('mm: Failed to create geometry',0, GEOM_ERR)
402
403      if(.not. geom_ncent(geom, ncent) )
404     &    call errquit("mm:geom_ncent",0,0)
405
406      if(.not.ma_push_get(mt_dbl,3*ncent,'ctmp',h_ctmp,i_ctmp))
407     & call errquit( pname//'Failed to allocate memory for ctmp',
408     & 3*ncent, MA_ERR)
409
410      if(.not. geom_cart_coords_get(geom,dbl_mb(i_ctmp)))
411     &    call errquit("mm:geom_cart_coords_get",0,0)
412
413
414      call util_convert_units("au","angstrom",scale)
415      call dscal(3*ncent, scale,dbl_mb(i_ctmp),1)
416      call mm_set_coord(nt,int_mb(i_iact),dbl_mb(i_ctmp))
417      call mm_vdw_qmcoords_load(rtdb)
418      call mm_vdw14_qmcoords_load(rtdb)
419      call mm_bond_qmcoords_load(rtdb)
420
421      if(.not.ma_pop_stack(h_ctmp))
422     & call errquit('mm:
423     &              Failed to deallocate stack c_tmp',ncent,
424     &       MA_ERR)
425
426      if(.not.geom_destroy(geom))
427     & call errquit('mm: Failed to destroy geometry',0, GEOM_ERR)
428
429
430      end
431
432      subroutine mm_set_coord(nt,ai,c)
433      implicit none
434#include "mafdecls.fh"
435#include "errquit.fh"
436#include "geom.fh"
437#include "rtdb.fh"
438#include "global.fh"
439#include "inp.fh"
440#include "mm_bond_coords_data.fh"
441#include "mm_link_data.fh"
442#include "mm_geom_data.fh"
443#include "mm_coords_data.fh"
444
445      integer rtdb
446      integer nt
447      integer ai(nt)
448      double precision c(3,nt)
449
450      integer i,j,k
451      integer indx, iqm, imm
452      character*32 pname
453
454      pname = "mm_set_coord"
455c     write(*,*) 'in ', pname
456
457      do i=1,nt
458        do j=1,nqml
459          indx = int_mb(i_iqml+j-1)
460          if(ai(i).eq.indx) then
461            do k = 1,3
462              dbl_mb(i_rqml+3*(j-1)+k-1) = c(k,i)
463            end do
464          end if
465        end do
466      end do
467
468c     -------------------
469c     update i_rqm array
470c     -------------------
471
472      do i=1,nqm
473        iqm = int_mb(i_iqm+i-1)
474        do j=1,nqml
475          indx = int_mb(i_iqml+j-1)
476          if(iqm.eq.indx) then
477            do k = 1,3
478              dbl_mb(i_rqm+3*(i-1)+k-1) = dbl_mb(i_rqml+3*(j-1)+k-1)
479            end do
480            exit
481          end if
482        end do
483      end do
484
485c     -------------------
486c     update i_rmm array
487c     -------------------
488      if(qmlink) then
489        do i=1,nlink
490          indx = int_mb(i_lb+nlink+i-1)
491          do j=1,nmm
492            imm = int_mb(i_imm+j-1)
493            if(imm.eq.indx) then
494              do k = 1,3
495                dbl_mb(i_rmm+3*(j-1)+k-1) = dbl_mb(i_rqml+
496     >                                      (nqm+i-1)*3+k-1)
497              end do
498              exit
499            end if
500          end do
501        end do
502      end if
503
504      end
505
506      subroutine mm_geom_restore(rtdb)
507      implicit none
508#include "mafdecls.fh"
509#include "errquit.fh"
510#include "geom.fh"
511#include "rtdb.fh"
512#include "mm_geom_data.fh"
513      integer rtdb
514      logical ignore
515      character*32 pname
516
517      pname = "mm_restore_geom"
518c     write(*,*) pname
519
520      if(.not.aux_geom) return
521
522      ignore = rtdb_delete(rtdb,'geometry')
523
524      if(oldgeom.ne.' ') then
525        if(.not.rtdb_cput(rtdb,'geometry',1,oldgeom))
526     &    call errquit(pname//' storing geom to rtdb', 0, RTDB_ERR)
527      end if
528
529      end
530
531      subroutine mm_geom_create_xyz_file(rtdb)
532      implicit none
533#include "util.fh"
534#include "mafdecls.fh"
535#include "errquit.fh"
536#include "rtdb.fh"
537#include "msgids.fh"
538#include "global.fh"
539#include "stdio.fh"
540#include "geom.fh"
541
542      integer rtdb
543
544      integer geom
545      character*30 pname
546      character*30 geomname
547      character*30 operation
548      character*50 filename
549      character*255 full_filename
550      integer ncent
551      logical status
552      logical ignore
553
554      pname = "mm_geom_create_xyz_file"
555      geomname = "geometry"
556
557
558      if (ga_nodeid().eq.0) then
559        call util_file_prefix('opt.xyz',filename)
560        call util_file_name_noprefix(filename,.false.,
561     >                                .false.,
562     >                                full_filename)
563
564        if(.not.geom_create(geom,'geom_tmp'))
565     >    call errquit('mm: Failed to create geometry',0, GEOM_ERR)
566        if(.not.geom_rtdb_load(rtdb,geom,geomname))
567     >    call errquit('mm: Failed to load geometry',0, GEOM_ERR)
568        status =  geom_ncent(geom,ncent)
569        open(88,file=full_filename,form='formatted')
570        if (.not. geom_print_xyz(geom, 88))
571     $    call errquit('mm:geom_print_xyz?',0, GEOM_ERR)
572        close(88)
573
574        if(.not.geom_destroy(geom))
575     >    call errquit('mm: Failed to destroy geometry',0, GEOM_ERR)
576
577      endif
578
579      end
580
581      subroutine mm_map_fixed_atoms(rtdb)
582      implicit none
583#include "rtdb.fh"
584#include "util.fh"
585#include "inp.fh"
586#include "mafdecls.fh"
587#include "errquit.fh"
588#include "geom.fh"
589#include "mm_link_data.fh"
590#include "mm_geom_data.fh"
591
592      integer rtdb
593
594      character*30 pname
595      integer i,j,k
596      integer nact0, nact1
597      integer i_cons,h_cons
598      integer i_act,h_act
599      integer ma_type
600      integer jj0,jj,ii,ii0
601      integer nal,h_al,i_al
602      integer h_am,i_am
603      logical aflag
604
605      pname = "mm_map_fixed_atoms"
606c     write(*,*) pname
607
608c     ---------------
609c     get active list
610c     ---------------
611      if (rtdb_ma_get(rtdb, 'geometry:actlist', ma_type,
612     $        nact0, h_cons)) then
613        if (.not. ma_get_index(h_cons, i_cons))
614     $           call errquit(pname,h_cons,
615     &       MA_ERR)
616
617        if (.not.rtdb_delete(rtdb, 'geometry:actlist'))
618     $       call errquit(pname,0,
619     &       RTDB_ERR)
620        do i=1,nact0
621        end do
622        if (.not.rtdb_put(rtdb, 'qmmm:actlist',
623     >       mt_int,nact0,int_mb(i_cons)))
624     $       call errquit(pname,0,
625     &       RTDB_ERR)
626
627      else
628        nact0 = nact
629        if(.not.ma_alloc_get(mt_int, nact0, 'qmmm actlist',
630     &      h_cons, i_cons) ) call errquit(
631     &      'qmmm_data_alloc: unable to allocate heap space',
632     &      nact, MA_ERR)
633        do i=1,nact0
634           int_mb(i_cons+i-1) = i
635        end do
636      end if
637
638c    ----------------------
639c    create new active list
640c    ----------------------
641      if(.not.ma_push_get(MT_INT, nfg, 'qmmm fixed atom list',
642     &      h_act, i_act) ) call errquit(
643     &      pname//' unable to allocate stack',
644     &      nact, MA_ERR)
645
646c     if(.not.ma_push_get(MT_LOG, ng, 'tmp qmmm act atom list',
647c    &      h_am, i_am) ) call errquit(
648c    &      pname//' unable to allocate stack',
649c    &      ng, MA_ERR)
650
651c     TP: call following subroutine to map active atom.
652c     For now, assume all (qm) atoms are active.
653
654c       call mm_cons_get_map(ng,int_mb(i_ig),
655c    >                       log_mb(i_am))
656
657c     --------------------------------------------
658c     find total number (nal) of link atoms
659c     that are bonded to active qm atoms
660c     and store their global index in i_al
661c     Note that link atoms are always indexed
662c     by the global index of the classical atom
663c     --------------------------------------------
664      if(nlink.ne.0) then
665
666        if(.not.ma_push_get(MT_INT, nlink, 'qmmm tmp link',
667     &        h_al, i_al) ) call errquit(
668     &        pname//' unable to allocate stack',
669     &        nact, MA_ERR)
670
671        nal=0
672        do i=1,nact0
673c         get global index next active atom in the auxilary geometry
674          ii0 = int_mb(i_cons+i-1)
675          ii  = int_mb(i_iqml+ii0-1)
676          aflag = .true.
677          if(aflag) then
678            do j=1,nlink
679c             get index of qm atom bonded to a link atom
680              jj  = int_mb(i_lb+j-1)
681c             if qm atom bonded to a link atom is active
682c             store corresponding index of the link atom
683              if(aflag.and.(jj.eq.ii)) then
684                nal=nal+1
685                int_mb(i_al+nal-1)=int_mb(i_lb+nlink+j-1)
686              end if
687            end do
688          end if
689        end do
690      else
691        nal = 0
692      end if
693
694c     ------------------------------------------------
695
696c     ------------------------------------------------
697c     construct active atom index
698c     note that if classical boundary atom is
699c     active the corresponding link atom is automatically
700c     flagged as active in the first do loop because it
701c     carries the same global index. The second do loop takes
702c     care of the case when qm boundary atom is active using
703c     link index constructed above
704c     ------------------------------------------------
705      nact1=0
706      do i=1,nfg
707        ii = int_mb(i_iqml+i-1)
708        do j=1,nact0
709         jj0 = int_mb(i_cons+j-1)
710         jj  = int_mb(i_iqml+jj0-1)
711         aflag = .true.
712         if(aflag.and.(ii.eq.jj)) then
713            nact1 = nact1 + 1
714            int_mb(i_act+nact1-1) = i
715            goto 1
716          end if
717        end do
718        do j=1,nal
719         jj  = int_mb(i_al+j-1)
720         if(ii.eq.jj) then
721            nact1 = nact1 + 1
722            int_mb(i_act+nact1-1) = i
723            goto 1
724          end if
725        end do
7261       continue
727      end do
728
729      if(nact1.ne.0) then
730        if (.not.rtdb_put(rtdb, 'geometry:actlist',
731     >          mt_int,nact1,int_mb(i_act)))
732     >         call errquit(pname,0,
733     >         RTDB_ERR)
734      end if
735
736      if (.not.ma_free_heap(h_cons))
737     $   call errquit(pname,h_cons,
738     &       MA_ERR)
739
740
741      if (.not. ma_chop_stack(h_act) ) call errquit(
742     &    pname//' ma_pop_stack ',
743     &    0, MA_ERR)
744
745      end
746
747      subroutine mm_restore_fixed_atoms(rtdb)
748      implicit none
749#include "rtdb.fh"
750#include "util.fh"
751#include "inp.fh"
752#include "mafdecls.fh"
753#include "errquit.fh"
754#include "nwc_const.fh"
755#include "geom.fh"
756      integer rtdb
757
758      character*30 pname
759      integer nact
760      integer i_cons,h_cons
761      integer ma_type
762      logical ignore
763
764      pname = "mm_restore_fixed_atoms"
765c     write(*,*) pname
766
767      ignore = rtdb_delete(rtdb, 'geometry:actlist')
768
769      if (rtdb_ma_get(rtdb, 'qmmm:actlist', ma_type,
770     $        nact, h_cons)) then
771            if (.not. ma_get_index(h_cons, i_cons))
772     $           call errquit(pname,h_cons,
773     &       MA_ERR)
774      else
775        return
776      end if
777
778      if (.not.rtdb_delete(rtdb, 'qmmm:actlist'))
779     $       call errquit(pname,0,
780     &       RTDB_ERR)
781
782      if (.not.rtdb_put(rtdb, 'geometry:actlist',
783     >       mt_int,nact,int_mb(i_cons)))
784     $       call errquit(pname,0,
785     &       RTDB_ERR)
786
787      if (.not.ma_free_heap(h_cons))
788     $   call errquit(pname,h_cons,
789     &       MA_ERR)
790
791      end
792