| 查看: 169 | 回复: 0 | |||
| 当前主题已经存档。 | |||
[交流]
【求助】计算库伦相互作用能的fortran程序并行化(1)
|
|||
|
c 此程序是化学方面理论计算部分用于计算静电相互作用能（即库仑作用能）和力的一个被调用子程序， 在计算能量e时为四个大的循环部分（这四部分可细分为9个小部分分别是循环的）： 原子部分（myid==1） 原子-原子 Atom-atom 原子-键 Atom-bond 原子-孤对电子 Atom-lp 原子-∏电子 Atom-pie 键部分（myid==2） 键-键 Bond-bond 键-∏电子 Bond lp 孤对电子部分（myid==3） 孤对电子-孤对电子 Lp-lp 孤对电子-∏电子 Lp-pie ∏电子部分（myid==4） ∏电子-∏电子 Pie-pie 现在我的想法就是把这四部分实现并行，请教大家。 （在我想分的进程处加了我的分配语句，已经注销了，请大家帮忙实现并行，有其他的办法更好） ############################################################### c ## ## c ## subroutine echarge1a -- charge derivs via double loop c ## ## c ############################################################### c c c "echarge1a" calculates the charge-charge interaction energy c and first derivatives with respect to Cartesian coordinates c using a pairwise double loop c c subroutine echarge1a implicit none c-----------------------------09.04.21 lc-------------------- include 'mpi.i' c-----------------------------09.04.21 lc-------------------- include 'sizes.i' include 'atoms.i' include 'bound.i' include 'cell.i' include 'cutoff.i' include 'cutoffa.i' include 'charge.i' include 'chgpot.i' include 'couple.i' include 'deriv.i' include 'energi.i' include 'group.i' include 'inter.i' include 'molcul.i' include 'shunt.i' include 'units.i' include 'usage.i' include 'virial.i' integer i,j,k,iii,jjj integer ii,kk integer in,kn integer ic,kc real*8 e,de,dc real*8 f,fi,fik real*8 r,r2,fgrp real*8 xi,yi,zi real*8 xr,yr,zr real*8 xc,yc,zc real*8 xic,yic,zic real*8 dedx,dedy,dedz real*8 dedxc,dedyc,dedzc real*8 shift,taper,trans real*8 dtaper,dtrans real*8 rc,rc2,rc3,rc4 real*8 rc5,rc6,rc7 real*8 vxx,vyy,vzz real*8 vyx,vzx,vzy real*8 cscale(maxatm) logical proceed,iuse real*8 dipolex,dipoley,dipolez integer at,istep,t,pi,M1,M2,MM1,MM2 real*8 xcl,ycl,zcl,xrl,yrl,zrl,r2l,kl,einterobj include 'atmtyp.i' include 'bath.i' include 'abeem.i' c-----------------------------09.04.21 zj-------------------- real*8 cputime11,cputime12,cputime21,cputime22 real*8 cputime31,cputime32,cputime41,cputime42 c real*8 cputime1,cputime2,cputime1,cputime2 c-----------------------------09.04.21 zj-------------------- c c c zero out the charge interaction energy and derivatives c c open(11,file='each-ec.txt',status='unknown') c open(12,file='ec.txt',status='unknown') c-----------------------------09.04.21 lc-------------------- open(123,file='time.txt',status='unknown') call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr) c-----------------------------09.04.21 lc-------------------- nion=n natom=nion nmol1=nmol c*******trans the coordinates into simutanious ones do i=1,natom ca(i,1)=x(i) ca(i,2)=y(i) ca(i,3)=z(i) an(i)=atomic(i) molcule1(i)=molcule(i) am(i)=an(i) nt12(i)=n12(i) do j=1,nt12(i) t12(j,i)=i12(j,i) enddo enddo call peptidecouple do i=1,na pchg(i)=bcharge(i) molcule(I)=molcule1(i) x(I)=ca(i,1) Y(I)=ca(i,2) z(i)=ca(i,3) C write(*,*) pchg(i),NATOM,NLP,NATOMBOND enddo c open(UNIT=8,file='dipole',status='unknown') c dipolex=0.0d0 c dipoley=0.0d0 c dipolez=0.0d0 c do i=1,na c if(molcule(I).eq.obj)then c dipolex=dipolex+ca(i,1)*pchg(i) c dipoley=dipoley+ca(i,2)*pchg(i) c dipolez=dipolez+ca(i,3)*pchg(i) c endif c enddo c dipole=DSQRT(dipolex**2+dipoley**2+dipolez**2) c &*2.541765/0.529 c write(8,*) 'dipole= ',dipole c write(*,*) dipole,'dipole' c zero out the charge interaction energy and derivatives c ec = 0.0d0 ec1=0.0d0 e=0.0d0 do i = 1, n dec(1,i) = 0.0d0 dec(2,i) = 0.0d0 dec(3,i) = 0.0d0 end do if (nion .eq. 0) return c c set conversion factor and switching function coefficients c f = electric / dielec f=332.05382d0 use_group=.false. call switch ('CHARGE') cut = chgcuta off =chgoffa c c compute the charge interaction energy and first derivatives c do ii=1,na iion(ii)=ii jion(ii)=ii kion(ii)=ii enddo c--------------------------09.04.21 ------------------------------ c if(myid==1.or.myid==2.or.myid==3.or.myid==4.or.myid==5)then c write(123,*)"process",myid,"is alive" c--------------------------09.04.21 ------------------------------ do ii = 1, nion c2scale= 0.0d0 c4scale= 0.57d0 c3scale= 0.0d0 i = iion(ii) in = jion(ii) ic = kion(ii) xic = x(ic) yic = y(ic) zic = z(ic) xi = x(i) - xic yi = y(i) - yic zi = z(i) - zic fi = f * pchg(ii) iuse = (use(i) .or. use(ic)) do j = 1, na if(ii.ne.j) cscale(j) = 0.57d0 end do DO J=1,NATOMLP(II) CSCALE(ATOMLP(J,II))=0.0D0 ENDDO DO J=1,NATOMPIE(II) CSCALE(ATOMPIE(J,II))=0.0D0 ENDDO do j = 1, n12(ii) cscale(i12(j,ii)) = c2scale if(atomlp(1,i12(j,in)).ne.0) cscale(atomlp(1,i12(j,in)))=c2scale if(atomlp(2,i12(j,in)).ne.0) cscale(atomlp(2,i12(j,in)))=c2scale if(atompie(1,i12(j,in)).ne.0)then cscale(atompie(1,i12(j,in)))=c2scale endif if(atompie(2,i12(j,in)).ne.0) then cscale(atompie(2,i12(j,in)))=c2scale endif end do do j = 1, natbond12(ii) cscale(atbond12(j,ii)) = c2scale end do do j = 1, n13(ii) cscale(i13(j,ii)) = c3scale if(atomlp(1,i13(j,in)).ne.0)then cscale(atomlp(1,i13(j,in)))=c3scale endiF if(atomlp(2,i13(j,in)).ne.0) cscale(atomlp(2,i13(j,in)))=c3scale if(atompie(1,i13(j,in)).ne.0) then cscale(atompie(1,i13(j,in)))=c3scale endif if(atompie(2,i13(j,in)).ne.0) then cscale(atompie(2,i13(j,in)))=c3scale endif end do do j = 1, natbond13(ii) cscale(atbond13(j,ii)) = c2scale end do do j = 1, n14(ii) cscale(i14(j,ii)) = c4scale end do do j = 1, natbond14(ii) cscale(atbond14(j,ii)) = c3scale end do c c decide whether to compute the current interaction c c-----------------------------09.04.21-------------------- c if(myid==1)then c------------------------------------------------------------ do kk = ii+1, nion k = iion(kk) kn = jion(kk) kc = kion(kk) proceed = .true. c if (proceed) proceed = (iuse .or. use(k) .or. use(kc)) c c compute the energy contribution for this interaction c if (proceed) then xc = xic - x(kk) yc = yic - y(kk) zc = zic - z(kk) if (use_image) call image (xc,yc,zc,0) rc2 = xc*xc + yc*yc + zc*zc xr = xc + xi - x(kk) + x(kk) yr = yc + yi - y(kk) + y(kk) zr = zc + zi - z(kk) + z(kk) r2 = xr*xr + yr*yr + zr*zr r = sqrt(r2) fik = fi * pchg(kk) * cscale(kk) if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then fik = fik*lamp2 endif e = fik / r c de = -fik / r2 dc = 0.0d0 if (molcule(ii) .ne. molcule(kk)) then einter = einter + e end if if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. & (molcule(ii).ne.obj.and.molcule(kk).eq.obj))then eintero= eintero + e end if if ((molcule(ii).eq.obj.and.molcule(kk).eq.obj))then eintrao= eintrao + e end if if(abs(e).gt.800)then write(*,*) 'echarge',ii,kk,molcule(II),molcule(kk) write(*,*) pchg(ii),pchg(kk),r endif c write(11,*) 'atom-atom',e ec1 = ec1 + e end if enddo c------------------------------------------------- do kk = natom+1, natombond fi = f * pchg(ii) k = iion(kk) kn = jion(kk) kc = kion(kk) proceed = .true. if (proceed) then xc = xic - x(kk) yc = yic - y(kk) zc = zic - z(kk) if (use_image) call image (xc,yc,zc,0) rc2 = xc*xc + yc*yc + zc*zc xr = xc + xi - x(kk) + x(kk) yr = yc + yi - y(kk) + y(kk) zr = zc + zi - z(kk) + z(kk) r2 = xr*xr + yr*yr + zr*zr r = sqrt(r2) fik = fi * pchg(kk) * cscale(kk) if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then fik = fik*lamp2 endif e = fik / r if (molcule(ii) .ne. molcule(bdatom12(1,kk))) then einter = einter + e end if if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then eintero = eintero + e end if if ((molcule(ii).eq.obj.and.molcule(kk).eq.obj))then eintrao = eintrao + e end if if(abs(e).gt.800)then write(*,*) ii,kk,molcule(II),molcule(kk) endif c write(11,*) 'atom-bond',e ec1=ec1+e end if end do c----------------------------------------------- do kk = natombond+1, NATOMBOND+NLP fi = f * pchg(ii) k = iion(kk) kn = jion(kk) kc = kion(kk) proceed = .true. xc = xic - x(kk) yc = yic - y(kk) zc = zic - z(kk) if (use_image) call image (xc,yc,zc,0) rc2 = xc*xc + yc*yc + zc*zc xr = xc + xi - x(kk) + x(kk) yr = yc + yi - y(kk) + y(kk) zr = zc + zi - z(kk) + z(kk) r2 = xr*xr + yr*yr + zr*zr r = sqrt(r2) if(cscale(kk).ne.0.0d0)then IF(AM1(II).GT.AM1(KK))THEN M1=AM1(II) M2=AM1(KK) ELSE M1=AM1(kk) M2=AM1(II) ENDIF IF(LPATOM(II).EQ.0.AND.LPATOM(KK).NE.0)THEN MM1=KK MM2=II ENDIF IF(LPATOM(II).NE.0.AND.LPATOM(KK).EQ.0)THEN MM1=II MM2=KK ENDIF if((M1.eq.182.and.M2.EQ.11).AND.distance(kk,ii).lt.2.0d0)then if(lamp2.ne.0.0d0)then cscale(kk)=0.6060d0-0.0496/(1.0d0+exp((r-1.3978)/0.1034)) endif endif if((M1.eq.881.and.M2.EQ.182).AND.distance(KK,iI).lt.2.0d0)then if(lamp2.ne.0.0d0)then cscale(KK)=0.6300d0-0.1297/(1.0d0+exp((r-1.8707)/0.2842)) endif endif if((M1.eq.883.and.M2.EQ.11).AND.distance(KK,iI).lt.2.0d0)then if(lamp2.ne.0.0d0)then cscale(KK)=0.63d0-0.0684/(1.0d0+exp((r-1.0099)/0.3376)) endif endif if((M1.eq.883.and.M2.EQ.881).AND.distance(kk,ii).lt.2.3d0)then if(lamp2.ne.0.0d0)then cscale(kk)=0.621d0-0.0813/(1.0d0+exp((r-1.151)/0.0697)) endif endif endif fik = fi * pchg(kk) * cscale(kk) c WRITE(*,*) AM(II),AM(KK),CSCALE(KK) if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then fik = fik*lamp2 endif e = fik / r if (molcule(ii) .ne. molcule(lpatom(kk))) then einter = einter + e end if if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then eintero = eintero + e end if if ((molcule(ii).eq.obj.and.molcule(kk).eq.obj))then eintrao = eintrao + e end if if(abs(e).gt.800)then write(*,*) ii,kk,molcule(II),molcule(kk) endif c write(11,*) 'atom-lp',e ec1=ec1+e end do c-------------------------------------------------- do kk = naTOMBOND+NLP+1,NA fi = f * pchg(ii) k = iion(kk) kn = jion(kk) kc = kion(kk) proceed = .true. xc = xic - x(kk) yc = yic - y(kk) zc = zic - z(kk) if (use_image) call image (xc,yc,zc,0) rc2 = xc*xc + yc*yc + zc*zc xr = xc + xi - x(kk) + x(kk) yr = yc + yi - y(kk) + y(kk) zr = zc + zi - z(kk) + z(kk) r2 = xr*xr + yr*yr + zr*zr r = sqrt(r2) fik = fi * pchg(kk) * cscale(kk) c WRITE(*,*) AM(II),AM(KK),CSCALE(KK) if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then fik = fik*lamp2 endif e = fik / r if (molcule(ii) .ne. molcule(PIEatom(kk))) then einter = einter + e end if if ((molcule(ii).eq.obj.and.molcule(kk).ne.obj).or. &(molcule(ii).ne.obj.and.molcule(kk).eq.obj))then eintero = eintero + e end if if ((molcule(ii).eq.obj.and.molcule(kk).eq.obj))then eintrao = eintrao + e end if if(abs(e).gt.800)then write(*,*) ii,kk,molcule(II),molcule(kk) endif c write(11,*) 'atom-pie',e ec1=ec1+e end do end do |
回复此楼» 猜你喜欢
寻求一种能扛住强氧化性腐蚀性的容器密封件
已经有3人回复
-
请问哪里可以有青B申请的本子可以借鉴一下。
已经有4人回复
-
真诚求助：手里的省社科项目结项要求主持人一篇中文核心，有什么渠道能发核心吗
已经有6人回复
-
孩子确诊有中度注意力缺陷
已经有14人回复
-
三甲基碘化亚砜的氧化反应
已经有4人回复
-
请问下大家为什么这个铃木偶联几乎不反应呢
已经有5人回复
-
请问有评职称，把科研教学业绩算分排序的高校吗
已经有5人回复
-
2025冷门绝学什么时候出结果
已经有3人回复
-
天津工业大学郑柳春团队欢迎化学化工、高分子化学或有机合成方向的博士生和硕士生加入
已经有4人回复
-
康复大学泰山学者周祺惠团队招收博士研究生
已经有6人回复
