24小时热门版块排行榜    

查看: 318  |  回复: 0

lxz73

铜虫 (小有名气)

[求助] 高手帮帮忙啊

这个程序应该是计算密度分布的 但是看不懂,请教高手
而且编译的时候还有警告!
c program calculate density profile

c  density profiles : (1) polymer; (2) nano; (3) solvent; (4) total
        parameter (NdpdMax=1400000)
        parameter (Ndenpro=500)
        parameter (nGn=15)

         implicit double precision (a-h,o-z)
        character*25 file_in,file_out

       real*8 rdpd(3,NdpdMax),vdpd(3,NdpdMax),rdpdcont(3,NdpdMax)
         real*8 denpro(4,Ndenpro)
         real*8 dz
         real*8 s_size(3),cell(3),sd,rho,D0
        Integer*4 species(NdpdMax)

        Integer*4 n_dpd_polymer,n_dpd_ps,n_dpd_nano
        integer*4 mdenpro
        Real*8 Ene(-1:3,-1:3),elsc,d_els,rm,T_solv,r_eq

                open(2, file='density.in')
              read(2,*)file_in, file_out, dz
            close(2)

            open(40, file=file_in)
                read(40,*)itt,ndpd,npoly,lchain,lchainA,n_dpd_ps
     &                ,n_dpd_nano
            read(40,*)T_solv
            read(40,*)ENE(-1,-1),ENE(1,1),ENE(-1,1),ENE(0,-1)
     &        , ENE(0,1),ENE(0,0),ENE(2,2),ENE(1,2)
          read(40,*)ENE(-1,2), ENE(0,2),ENE(3,3),ENE(2,3)
     &        ,  ENE(1,3),ENE(-1,3),ENE(0,3)
            read(40,*)s_size,D0

                D0=D0+2.75
                dv=s_size(1)*s_size(2)*dz

                n_dpd_polymer=npoly*lchain

                mdenpro=D0/dz+0.5
            do jn=1,mdenpro
                        denpro(1,jn)=0.0
                        denpro(2,jn)=0.0
                        denpro(3,jn)=0.0
                        denpro(4,jn)=0.0
            enddo
         
                do i=1,n_dpd_ps
            read(40,*)rdpd(1,i),rdpd(2,i),rdpd(3,i)
     >           ,vdpd(1,i),vdpd(2,i),vdpd(3,i)
     &           ,rdpdcont(1,i),rdpdcont(2,i),rdpdcont(3,i)
     &           ,species(i)

           m=rdpd(3,i)/dz+1

          denpro(4,m)=denpro(4,m)+1


         if(i.le.n_dpd_polymer) denpro(1,m)=denpro(1,m)+1
         if(i.gt.n_dpd_polymer.and.i.le.(n_dpd_polymer+n_dpd_nano)) then
           denpro(2,m)=denpro(2,m)+1
         endif
         if(i.gt.(n_dpd_polymer+n_dpd_nano).and.i.le.n_dpd_ps) then
           denpro(3,m)=denpro(3,m)+1
         endif

         enddo
                       
           open(13, file=file_out)
        do k=1,mdenpro
             write(13, '(f10.3,4f18.7)')(k-0.5)*dz,
     &        denpro(1,k)/dv,
     &            denpro(2,k)/dv,
     &        denpro(3,k)/dv,
     &        denpro(4,k)/dv
        enddo



       stop
        end
回复此楼

» 猜你喜欢

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子,希望有用哦~

科研从小木虫开始,人人为我,我为人人
相关版块跳转 我要订阅楼主 lxz73 的主题更新
信息提示
请填处理意见