1 2 subroutine smd_exlist_init_system() 3 implicit none 4#include "errquit.fh" 5#include "inp.fh" 6#include "mafdecls.fh" 7#include "rtdb.fh" 8#include "util.fh" 9#include "global.fh" 10c 11 character*32 sp_exlist,sp_atom 12 character*32 tag,pname 13 logical result 14 15 pname = "smd_exlist_init_system" 16c 17 tag = "atom" 18 call smd_system_get_component(sp_atom,tag,result) 19 if(.not.result) 20 > call errquit( 21 > pname//'no component '//tag,0,0) 22 23 tag = "excl_list" 24 call smd_system_get_component(sp_exlist,tag,result) 25 if(.not.result) 26 > call errquit( 27 > pname//'no component '//tag,0,0) 28 29 30 call smd_exlist_init(sp_exlist,sp_atom,result) 31c 32 if(.not.result) then 33 tag = "excl_list" 34 call smd_system_unset_component(tag) 35 end if 36 37 return 38 end 39 40 subroutine smd_exlist_init(sp_exlist,sp_atom,result) 41 implicit none 42#include "errquit.fh" 43#include "inp.fh" 44#include "mafdecls.fh" 45#include "rtdb.fh" 46#include "util.fh" 47#include "global.fh" 48c 49 character*(*) sp_exlist 50 character*(*) sp_atom 51 logical result 52c 53 character*32 pname 54 character*80 tag 55 integer np,nl 56 integer i_p,i_l,i_ir 57 integer h_l 58 integer i_list 59 integer i 60c 61 pname = "smd_exlist_init" 62c 63c write(*,*) "in "//pname 64c 65c get total number of atoms 66c ------------------------- 67 call smd_atom_ntot(sp_atom,np) 68c 69c gestimate the size of pair list 70c ------------------------------ 71 nl = min( np*200, ma_inquire_avail(MT_INT)) 72c 73c create pointer memory 74c --------------------- 75 call smd_namespace_create(sp_exlist) 76 tag = "exlist:pointer" 77 call smd_data_create(sp_exlist,"exlist:pointer",np,MT_INT) 78 call smd_data_get_index(sp_exlist,tag,i_p,result) 79 if(.not. result) 80 > call errquit( 81 > pname//'error getting index for'//tag,0, RTDB_ERR) 82 83c 84c create temporary scratch array for list since 85c we do not know the size yet 86c --------------------------------------------- 87 if(.not.ma_push_get(mt_int,nl,'tmp l',h_l,i_l)) 88 + call errquit(pname//'Failed to allocate memory for tmp l', 89 + nl, MA_ERR) 90 91c 92c get atom residue index for exclusion criteria 93c --------------------------------------------- 94 tag = "atom:resid" 95 call smd_data_get_index(sp_atom,tag,i_ir,result) 96 if(.not. result) 97 > call errquit( 98 > pname//'error getting index for'//tag,0, RTDB_ERR) 99 100 101 call smd_exlist_set(nl, 102 + int_mb(i_l), 103 + np, 104 + int_mb(i_p), 105 + int_mb(i_ir)) 106 107c 108c create list memory 109c nl now contains the actual size 110c ------------------------------- 111 112 if(nl.eq.0) then 113 result = .false. 114 call smd_namespace_destroy(sp_exlist) 115 goto 200 116 end if 117c 118c write(*,*) "excluded size list",nl 119 tag = "exlist:list" 120 call smd_data_create(sp_exlist,tag,nl,MT_INT) 121 call smd_data_get_index(sp_exlist,tag,i_list,result) 122 if(.not. result) 123 > call errquit( 124 > pname//'error getting index for'//tag,0, RTDB_ERR) 125 126 127 do i=1,nl 128 int_mb(i_list+i-1) = int_mb(i_l+i-1) 129 end do 130 131 if(.not.ma_pop_stack(h_l)) 132 & call errquit(pname//'Failed to deallocate stack h_l',nl, 133 & MA_ERR) 134 135200 continue 136 return 137 end 138c 139c 140 subroutine smd_exlist_set(nl,list,natms,point,resid) 141c 142c constructs excluded pairt list 143c point(i) is an index to list() array 144c that contains all the pairs of atom i 145c In other words all the atoms paired with atom i 146c are contained in list(point(i)),...,list(point(i+1)-1) 147c 148 implicit none 149#include "errquit.fh" 150 integer nl,natms 151 integer list(nl) 152 integer point(natms) 153 integer resid(natms) 154c 155 integer i,j 156 integer nlist 157 158 character*30 pname 159 160 pname = "smd_exlist_set" 161 162 nlist = 0 163 do i=1,natms-1 164 165 point(i)=nlist+1 166 167 do j=i+1,natms 168 169 if(resid(i).eq.resid(j))then 170 171 nlist=nlist+1 172 173 if(nlist.gt.nl) 174 + call errquit(pname//'exceeded list dimensions', 175 + nl, 0) 176 177 list(nlist)=j 178 179 endif 180 181 enddo 182 enddo 183c 184c we need to set this to mark 185c the end of pair list index for natms-1 186c it would be used as point(natms)-1 187c ------------------------------------- 188 point(natms)=nlist+1 189 190 nl = nlist 191 192 return 193 194 END 195c $Id$ 196