1 subroutine argos_cafe_polar(lself,lpbcs, 2 + xw,xwm,fw,pw,pwp,zwi,rtos,iwdt,iwz, 3 + iwfr,iwto,jwfr,jwto,xs,xsm,fs,ps,psp,zs, 4 + isga,isat,isdt,isgr,ismf,isml,isss,isq1,isq2,isq3,isgm,ishop,isz, 5 + isfr,isto,jsfr,jsto,lpbc,lstptr,lupden,eww,esw,ess,fss,esa,lseq) 6c 7c in log : lself : true for box self interactions 8c in r*8 : xw(mwm,3,mwa) : solvent coordinates 9c in r*8 : xwm(mwm,3) : solvent molecule center of mass coordinates 10c in/out r*8 : fw(nwm,3,mwa,2) : solvent forces 11c out r*8 : rtos(mwm) : solvent distances to closes solute atom 12c in int : iwdt(mwm) : solvent dynamics type 13c out int : iwz(mwm) : solvent boundary type 14c in int : iwfr,iwto : first and last solvent molecule i 15c in int : jwfr,jwto : first and last solvent molecule j 16c in r*8 : xs(msa,3) : solute atom coordinates 17c in r*8 : xsm(msm,3) : solute molecule center of mass coordinates 18c in/out r*8 : fs(msa,3,2) : solute forces 19c out r*8 : zs(msm,3) : solute virial 20c in int : isga(msa) : solute global atom number 21c in int : isat(msa) : solute atom type 22c in int : isdt(msa) : solute dynamics type 23c in int : isgr(msa) : solute charge group 24c in int : ismf(msa) : solute molecule fraction 25c in int : isml(msa) : solute molecule 26c in int : isss(msa) : solute separation shifted scaling type 27c in int : isq1(msa) : solute charge type 1 28c in int : isq2(msa) : solute charge type 2 29c in int : isq3(msa) : solute charge type 3 30c out int : isz(msa) : solute boundary type 31c in int : isfr,isto : first and last solute atom i 32c in int : jsfr,jsto : first and last solute atom j 33c in log : lpbc : flag for periodic boundary conditions 34c in/out int : lstptr : list pointer 35c in/out r*8 : eww(mpe,2) : mpe=1 solvent bond energy 36c 2 solvent angle energy 37c 3 solvent torsion energy 38c 4 solvent improper energy 39c 5 solvent van der Waals energy 40c 6 solvent electrostatic energy 41c 7 solvent potential energy 42c eswl(msf,2) : solute-solvent van der Waals energy 43c in/out r*8 : eswq(msf,2) : solute-solvent electrostatic energy 44c in/out r*8 : essl(msf,msf,2) : solute-solute van der Waals energy 45c in/out r*8 : essq(msf,msf,2) : solute-solute electrostatic energy 46c in/out r*8 : essr(msf,msf,2) : solute-solute reaction field energy 47c in/out r*8 : esp(msf,2) : solute potential energy 48c in/out r*8 : esb(msf,2) : solute bond energy 49c in/out r*8 : esh(msf,2) : solute angle energy 50c in/out r*8 : esd(msf,2) : solute torsion energy 51c in/out r*8 : eso(msf,2) : solute improper energy 52c in log : lupden : update bonded energy contrbutions 53c 54c dimensions nwm,nwa and nsa need to have been given by a call to argos_cafe_initx 55c 56 implicit none 57c 58#include "argos_cafe_common.fh" 59#include "mafdecls.fh" 60c 61 real*8 xw(mwm,3,mwa),xwm(mwm,3),fw(mwm,3,mwa,2),zwi(3,3,2) 62 real*8 xs(msa,3),xsm(msm,3),fs(msa,3,2) 63 real*8 pw(mwm,3,mwa,2),pwp(mwm,3,mwa,2,2) 64 real*8 ps(msa,3,2),psp(msa,3,2,2) 65 integer iwdt(mwm),iwz(mwm),isz(msa) 66 integer isga(msa),isat(msa),isdt(msa),isgr(msa),ismf(msa) 67 integer isml(msa),isss(msa),isq1(msa),isq2(msa),isq3(msa) 68 integer isgm(msa),ishop(msa) 69 integer lseq(mseq) 70 integer iwfr,iwto,jwfr,jwto,isfr,isto,jsfr,jsto 71 integer lstptr 72 logical lself,lpbc,lupden,lpbcs,lforce 73 real*8 rtos(mwm),zs(msf,3,32) 74 real*8 eww(mpe,2),esw(msf,mpe,2),ess(msf,msf,mpe,2),esa(nsa) 75 real*8 fss(msf,msf,3,2) 76c 77 integer nwloc,nsloc,nwnon,nsnon,npairs 78 integer lptr 79 integer i,j 80c 81c print*,'CAFE forces' 82c 83 do 1 j=1,2 84 do 2 i=1,3 85 zw(1,i,j)=zwi(1,i,j) 86 zw(2,i,j)=zwi(2,i,j) 87 zw(3,i,j)=zwi(3,i,j) 88 2 continue 89 1 continue 90c 91 if(lself) then 92 jwfr=iwfr 93 jwto=iwto 94 jsfr=isfr 95 jsto=isto 96 endif 97c 98 nwloc=iwto-iwfr+1 99 if(iwfr.eq.0.or.iwto.lt.iwfr) nwloc=0 100 nwnon=jwto-jwfr+1 101 if(jwfr.eq.0.or.jwto.lt.jwfr) nwnon=0 102 nsloc=isto-isfr+1 103 if(isfr.eq.0.or.isto.lt.isfr) nsloc=0 104 nsnon=jsto-jsfr+1 105 if(jsfr.eq.0.or.jsto.lt.jsfr) nsnon=0 106c 107c list [ list(1)=int_mb(lptr) ] with offset ndxp 108c 109c list( 1) : index to pairlist for solvent-solvent 110c list( 2) : index to pairlist for solute-solvent 111c list( 3) : index to pairlist for solvent-solute 112c list( 4) : index to pairlist for solute-solute 113c list( 5) : index to pairlist for solute bonds 114c list( 6) : index to pairlist for solute angles 115c list( 7) : index to pairlist for solute dihedrals 116c list( 8) : index to pairlist for solute impropers 117c list( 9) : index to pairlist for solute third neighbors 118c list(10) : index to pairlist for solute excluded pairs 119c 120c 121c pairlists 122c --------- 123c 124 if(lpair) call argos_pairs(lself,lpbcs,xw,xwm,iwdt,iwz, 125 + iwfr,iwto,jwfr,jwto,xs,xsm, 126 + isga,isat,isdt,isgr,isgm,ismf,isml,isss,isq1,isq2,isq3,ishop,isz, 127 + isfr,isto,jsfr,jsto,lpbc,lstptr,lseq) 128c 129 if(.not.lforce) return 130c 131c forces 132c ------ 133c 134 ndxp=lstptr 135 lptr=i_list+ndxp+24 136c 137 if(lself.and.nwloc.gt.0) then 138 call argos_cafe_fw(iwfr,iwto,xw,fw,iwdt,int_mb(i_iwa), 139 + int_mb(i_iwq), 140 + lpbc,eww,dbl_mb(i_vdw),dbl_mb(i_chg),mbt(1),numb(1),mbp(1), 141 + int_mb(i_ibnd(1)),dbl_mb(i_bnd(1)),dbl_mb(i_rbnd(1)), 142 + mht(1),numh(1),mhp(1), 143 + int_mb(i_iang(1)),dbl_mb(i_ang(1)),dbl_mb(i_rang(1)), 144 + dbl_mb(i_rub(1)), 145 + mdt(1),numd(1),mdp(1), 146 + int_mb(i_idih(1)),dbl_mb(i_dih(1)),dbl_mb(i_rdih(1)), 147 + mit(1),numi(1),mip(1), 148 + int_mb(i_iimp(1)),dbl_mb(i_imp(1)),dbl_mb(i_rimp(1)), 149 + mtt(1),numt(1),int_mb(i_itrd(1)), 150 + mxt(1),numx(1),int_mb(i_ixcl(1))) 151 endif 152c 153c solvent-solvent forces 154c 155 npairs=int_mb(lptr) 156 if(npairs.gt.0) then 157 call argos_cafe_fpww(xw,xwm,fw,pw,pwp,iwdt,iwfr,nwloc,lpbc,eww, 158 + dbl_mb(i_vdw),dbl_mb(i_chg),int_mb(i_iwa),int_mb(i_iwq), 159 + int_mb(i_s2i1), 160 + int_mb(lptr+1),int_mb(lptr+1+2*nwloc),int_mb(lptr+1+4*nwloc), 161 + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1), 162 + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4), 163 + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5), 164 + dbl_mb(i_smr5),dbl_mb(i_smr6)) 165 lptr=lptr+4*nwloc+1+npairs 166 else 167 lptr=lptr+1 168 endif 169c 170c solute-solvent forces 171c 172 npairs=int_mb(lptr) 173 if(npairs.gt.0) then 174 call argos_cafe_fpsw(xs,xsm,fs,zs,ps,psp,isga,isat,isdt, 175 + ismf,isml,isss, 176 + isq1,isfr,nsloc, 177 + xw,xwm,fw,pw,pwp,rtos,iwdt,lpbc,lpbcs,esw,esa, 178 + dbl_mb(i_vdw),dbl_mb(i_chg),int_mb(i_iwa),int_mb(i_iwq), 179 + int_mb(i_ias),int_mb(i_s2i1), 180 + int_mb(lptr+1),int_mb(lptr+1+2*nsloc),int_mb(lptr+1+4*nsloc), 181 + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1), 182 + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4), 183 + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),dbl_mb(i_s1r6), 184 + int_mb(i_s1i1),int_mb(i_s1i2),int_mb(i_s1i3), 185 + dbl_mb(i_smr5),dbl_mb(i_smr6)) 186 lptr=lptr+4*nsloc+1+npairs 187 else 188 lptr=lptr+1 189 endif 190c 191c solvent-solute forces 192c 193 npairs=int_mb(lptr) 194 if(npairs.gt.0) then 195 call argos_cafe_fpsw(xs,xsm,fs,zs,ps,psp,isga,isat, 196 + isdt,ismf,isml,isss, 197 + isq1,jsfr,nsnon, 198 + xw,xwm,fw,pw,pwp,rtos,iwdt,lpbc,lpbcs,esw,esa, 199 + dbl_mb(i_vdw),dbl_mb(i_chg),int_mb(i_iwa),int_mb(i_iwq), 200 + int_mb(i_ias),int_mb(i_s2i1), 201 + int_mb(lptr+1),int_mb(lptr+1+2*nsnon),int_mb(lptr+1+4*nsnon), 202 + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1), 203 + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4), 204 + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),dbl_mb(i_s1r6), 205 + int_mb(i_s1i1),int_mb(i_s1i2),int_mb(i_s1i3), 206 + dbl_mb(i_smr5),dbl_mb(i_smr6)) 207 lptr=lptr+4*nsnon+1+npairs 208 else 209 lptr=lptr+1 210 endif 211c 212c solute-solute forces 213c 214 npairs=int_mb(lptr) 215 if(npairs.gt.0) then 216 call argos_cafe_fpss(xs,xsm,fs,zs,ps,psp,isga,isat, 217 + isdt,ismf,isml,isss, 218 + isq2,isq3,isfr,nsloc,lpbc,lpbcs,ess,fss,esa, 219 + dbl_mb(i_vdw),dbl_mb(i_chg), 220 + int_mb(i_ias),int_mb(i_s2i1), 221 + int_mb(lptr+1),int_mb(lptr+1+2*nsloc),int_mb(lptr+1+4*nsloc), 222 + dbl_mb(i_smr1),dbl_mb(i_smr2),dbl_mb(i_s3r1),dbl_mb(i_s1r1), 223 + dbl_mb(i_s1r2),dbl_mb(i_s1r3),dbl_mb(i_s3r2),dbl_mb(i_s1r4), 224 + dbl_mb(i_smr3),dbl_mb(i_smr4),dbl_mb(i_s1r5),dbl_mb(i_s1r6), 225 + int_mb(i_s1i1),int_mb(i_s1i2),int_mb(i_s1i3),int_mb(i_s1i4), 226 + int_mb(i_s1i5),dbl_mb(i_smr5),dbl_mb(i_smr6), 227 + dbl_mb(i_smr5),dbl_mb(i_smr6)) 228 lptr=lptr+4*nsloc+1+npairs 229 else 230 lptr=lptr+1 231 endif 232c 233c solute bond forces 234c 235 npairs=int_mb(lptr) 236 if(npairs.gt.0) then 237 call argos_cafe_fsb(npairs,int_mb(lptr+1),mbt(2),mbp(2), 238 + int_mb(i_ibnd(2)),dbl_mb(i_bnd(2)),dbl_mb(i_rbnd(2)), 239 + max(isto,jsto),msa,isga,isgm,ismf,isdt,isq1,dbl_mb(i_chg), 240 + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq)) 241 lptr=lptr+1+npairs 242 else 243 lptr=lptr+1 244 endif 245c 246c solute angle forces 247c 248 npairs=int_mb(lptr) 249 if(npairs.gt.0) then 250 call argos_cafe_fsh(npairs,int_mb(lptr+1),mht(2),mhp(2), 251 + int_mb(i_iang(2)),dbl_mb(i_ang(2)),dbl_mb(i_rang(2)), 252 + dbl_mb(i_rub(2)), 253 + max(isto,jsto),msa,isga,isgm,ismf,isdt,isq1,dbl_mb(i_chg), 254 + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq)) 255 lptr=lptr+1+npairs 256 else 257 lptr=lptr+1 258 endif 259c 260c solute torsion forces 261c 262 npairs=int_mb(lptr) 263 if(npairs.gt.0) then 264 call argos_cafe_fsd(npairs,int_mb(lptr+1),mdt(2),mdp(2), 265 + int_mb(i_idih(2)),dbl_mb(i_dih(2)),dbl_mb(i_rdih(2)), 266 + max(isto,jsto),msa,isga,isgm,ismf,isdt, 267 + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq)) 268 lptr=lptr+1+npairs 269 else 270 lptr=lptr+1 271 endif 272c 273c solute improper dihedral forces 274c 275 npairs=int_mb(lptr) 276 if(npairs.gt.0) then 277 call argos_cafe_fso(npairs,int_mb(lptr+1),mit(2),mip(2), 278 + int_mb(i_iimp(2)),dbl_mb(i_imp(2)),dbl_mb(i_rimp(2)), 279 + max(isto,jsto),msa,isga,isgm,ismf,isdt, 280 + xs,fs,ess,lpbc,lpbcs,lupden,.true.,dbl_mb(i_sti),int_mb(i_lseq)) 281 lptr=lptr+1+npairs 282 else 283 lptr=lptr+1 284 endif 285c 286c solute third neighbor forces 287c 288 npairs=int_mb(lptr) 289 if(npairs.gt.0) then 290 call argos_cafe_fst(npairs,int_mb(lptr+1),mtt(2), 291 + int_mb(i_itrd(2)), 292 + dbl_mb(i_vdw),dbl_mb(i_chg), 293 + max(isto,jsto),msa,isat,isga,isgm,ismf,isdt,isq3,isss, 294 + xs,fs,ess,lpbc,lpbcs,dbl_mb(i_sti),esa,int_mb(i_lseq)) 295 lptr=lptr+1+npairs 296 else 297 lptr=lptr+1 298 endif 299c 300 ndxp=lptr-i_list 301c 302 do 3 j=1,2 303 do 4 i=1,3 304 zwi(1,i,j)=zw(1,i,j) 305 zwi(2,i,j)=zw(2,i,j) 306 zwi(3,i,j)=zw(3,i,j) 307 4 continue 308 3 continue 309c 310 return 311 end 312c $Id$ 313