24小时热门版块排行榜    

查看: 1607  |  回复: 7
本帖产生 1 个 程序强帖 ,点击这里进行查看

kathy2008

木虫 (正式写手)

[求助] 如何从高斯输出文件快速提出 pai 轨道信息。

如题。从高斯输出文件提出了eigenvector那一部分出来,即附件1。现在需要得到 pai 轨道信息。即附件2。 附件2 对应于附件1的32号,35号,38号,39号,40号,41号轨道(占据轨道),42号一直到47号(非占据轨道)的2Px值。求一小程序。请指点。谢谢。
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

【答案】应助回帖

★ ★
ben_ladeng(金币+2): 很详细,待楼主评定后奖励程序强帖 2011-06-18 17:42:21
kathy2008(金币+10): 2011-06-19 13:00:43
微尘、梦想(程序强帖+1): 2011-06-19 17:04:28
大概这个样子。只是需要手工输入轨道号(这样可能灵活一些),每次输入一个轨道序号,回车,输入 0 则结束整个程序……
CODE:
program ei
real, dimension(:,:), allocatable :: px,ppx
character(len=256) :: line
character(len=40) :: fm
integer :: nrow, ncol, i, j, k, ios

open(unit=12, file='eigenvector.out', status='old')
open(unit=13, file='2px.out', status='new')

do
   read(12,'(a)', iostat=ios) line
   if (ios /= 0) exit
   if (index(line,'EIGENVALUES') /= 0) then
      nrow=0
      ncol=0
      do
         read(12,'(a)', iostat=ios) line
         if (ios /= 0) exit
         if (line(1:4) == '    ') exit
         ncol=ncol+1
         if (index(line, '2PX') /= 0) nrow=nrow+1
      end do
      exit
   end if
end do

!write (*,*) nrow, ncol
rewind (12)

allocate(px(nrow,ncol),ppx(nrow,ncol))

i=0
j=0
do
   read(12,'(a)', iostat=ios) line
   if (ios /= 0) exit
   if (i == nrow) then
       i=0
       j=j+n
   end if
   if (index(line, '2PX') /= 0) then
       line = line(21:)
!      write(*,*) trim(line)
       i=i+1
       n = len_trim(line)/10
       write(fm,'(a,i0,a)') '(', n, 'f10.5)'
!      write(*,*) j
       read(line,fm) px(i,(j+1):(j+n))
   end if
end do

k=0
do
   write(*,*) 'please input a number between 1 and ', nrow, 'end the program by 0.'
   read(*,*) i
   if(i==0) exit
   k=k+1
   ppx(:,k) = px(:,i)
end do

!write(*,*) k/5, mod(k,5)

if (k>=5) then
   do j=1,k/5
      do i=1,nrow
         write(13,'(5f10.5)') ppx(i,(j-1)*5+1:j*5)
      end do
      write(13,*)
   end do
end if
if (mod(k,5) /=0) then
   write(fm,'(a,i0,a)') '(', mod(k,5), 'f10.5)'
   do i=1,nrow
      write(13, fm) ppx(i,(k/5*5+1):k)
   end do
end if

end program ei

2楼2011-06-18 16:41:08
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kathy2008

木虫 (正式写手)


dubo(金币+1): 欢迎常来程序语言版讨论 2011-07-31 13:36:05
引用回帖:
Originally posted by snoopyzhao at 2011-06-18 16:41:08:
大概这个样子。只是需要手工输入轨道号(这样可能灵活一些),每次输入一个轨道序号,回车,输入 0 则结束整个程序……

[code]
program ei
real, dimension(:,, allocatable :: px,ppx
character(len=256 ...

利用该程序提取π轨道信息,报错。信息如下
At line 48 of file eigen-pai-nc3h7-r2.f
Fortran runtime error: Bad value during floating point read

google也没有找到应对之策,请高手指点。
3楼2011-07-03 15:06:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)


dubo(金币+1): 欢迎常来程序语言版讨论 2011-07-31 13:36:26
你的 .out 文件是咋生成的,这次的这个文件比上次的文件每一行前面多了一个空格……

所以,你把程序中:
CODE:
line = line(21:)

改成
CODE:
line = line(22:)

就可以了……
4楼2011-07-03 22:28:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)


dubo(金币+1): 欢迎常来程序语言版讨论 2011-07-31 13:36:35
这样可能更好一些:
CODE:
      program ei
      implicit none
      real, dimension(:,:), allocatable :: px,ppx
      character(len=256) :: line
      character(len=40) :: fm
      integer :: nrow, ncol, i, j, k, ios, n, m
      
      open(unit=12, file='nc3h7-r2-sto-eiv.out', status='old')
      open(unit=13, file='nc3h7-r2-sto-pai.out', status='new')
      
      do
         read(12,'(a)', iostat=ios) line
         if (ios /= 0) exit
         if (index(line,'Eigenvalues') /= 0) then
            nrow=0
            ncol=0
            do
               read(12,'(a)', iostat=ios) line
               if (ios /= 0) exit
               if (line(1:4) == '    ') exit
               ncol=ncol+1
               if (index(line, '2PZ') /= 0) nrow=nrow+1
            end do
            exit
         end if
      end do
      
      write (*,*) nrow, ncol
      rewind (12)
      
      allocate(px(nrow,ncol),ppx(nrow,ncol))
      
      i=0
      j=0
      do
         read(12,'(a)', iostat=ios) line
         if (ios /= 0) exit
         if (i == nrow) then
             i=0
             j=j+n
         end if
         if (index(line, '2PZ') /= 0) then
             line = line(21:)
!            write(*,*) trim(line)
             i=i+1
             n = len_trim(line)/10
             m = mod(len_trim(line),10)
!            write (*,*) m, n
             if (m /= 0) then
                write(fm,'(a,i0,a,i0,a)') '(tr',m,',',n,'f10.5)'
             else
                write(fm,'(a,i0,a)') '(',n,'f10.5)'
             end if
!            write (*,*) fm
!            write(*,*) j
             read(line,fm) ppx(i,(j+1):(j+n))
         end if
      end do
      
      k=0
      do
         write(*,*) 'please input a number between 1 and ',nrow,',
     & end the program by 0.'
         read(*,*) i
         if(i==0) exit
         k=k+1
         ppx(:,k) = px(:,i)
      end do
      
      !write(*,*) k/5, mod(k,5)
      
      if (k>=5) then
         do j=1,k/5
            do i=1,nrow
               write(13,'(5f10.5)') ppx(i,(j-1)*5+1:j*5)
            end do
            write(13,*)
         end do
      end if
      if (mod(k,5) /=0) then
         write(fm,'(a,i0,a)') '(', mod(k,5), 'f10.5)'
         do i=1,nrow
            write(13, fm) ppx(i,(k/5*5+1):k)
         end do
      end if
      
      end program ei

5楼2011-07-03 22:59:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kathy2008

木虫 (正式写手)


dubo(金币+1): 欢迎常来程序语言版讨论 2011-07-31 13:36:56
引用回帖:
Originally posted by snoopyzhao at 2011-07-03 22:28:17:
你的 .out 文件是咋生成的,这次的这个文件比上次的文件每一行前面多了一个空格……

所以,你把程序中:
CODE:
line = line(21:)

改成
CODE:
line = line(22:)

就可以了……

运行后报错信息如下:
*** glibc detected *** ./a.out: double free or corruption (out): 0x00000000079e2980 ***
======= Backtrace: =========
/lib64/libc.so.6[0x3e64871ce2]
/lib64/libc.so.6(cfree+0x8c)[0x3e6487590c]
./a.out[0x4012c5]
./a.out[0x4019de]
/lib64/libc.so.6(__libc_start_main+0xf4)[0x3e6481d974]
./a.out[0x400a99]
6楼2011-07-04 07:23:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)


jjdg(金币+1): 感谢参与 2011-07-04 12:46:36
我才注意到,你把我在二楼给出的程序中的
CODE:
read(line,fm) px(i,(j+1):(j+n))

改成了
CODE:
read(line,fm) ppx(i,(j+1):(j+n))

害得我弄了半天才知道为啥出来的结果总是不对。
7楼2011-07-04 09:48:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)

★ ★
jjdg(金币+2): 辛苦了 2011-07-04 12:46:22
另外,原程序中
CODE:
write(*,*) 'please input a number between 1 and ', nrow, 'end the program by 0.'

应该改为
CODE:
write(*,*) 'please input a number between 1 and ', ncol, 'end the program by 0.'

8楼2011-07-04 09:50:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 kathy2008 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 华东师范大学-071000生物学-293分-求调剂 +3 研究生何瑶明 2026-03-18 3/150 2026-03-21 01:30 by JourneyLucky
[考研] 280求调剂 +7 咕噜晓晓 2026-03-18 8/400 2026-03-21 01:27 by JourneyLucky
[考研] 一志愿华南师大 070300(化学)304分求调剂 +3 0703武芊慧雪304 2026-03-18 3/150 2026-03-21 00:48 by JourneyLucky
[考研] 288求调剂 +16 于海海海海 2026-03-19 16/800 2026-03-20 22:28 by JourneyLucky
[考研] 290求调剂 +7 ^O^乜 2026-03-19 7/350 2026-03-20 21:43 by JourneyLucky
[考研] 一志愿武理材料工程348求调剂 +3  ̄^ ̄゜汗 2026-03-19 4/200 2026-03-20 21:01 by zhukairuo
[考研] 一志愿北京化工大学0703化学318分,有科研经历,求调剂 +4 一瓶苯甲酸 2026-03-14 4/200 2026-03-20 20:36 by fen_rao
[考研] 一志愿 南京航空航天大学大学 ,080500材料科学与工程学硕 +5 @taotao 2026-03-20 5/250 2026-03-20 20:16 by JourneyLucky
[考研] 广西大学家禽遗传育种课题组2026年硕士招生(接收计算机专业调剂) +3 123阿标 2026-03-17 3/150 2026-03-20 15:58 by 飞行琦
[考博] 招收博士1-2人 +3 QGZDSYS 2026-03-18 3/150 2026-03-20 11:58 by 呱呱呱呱叫
[考研] 材料学硕318求调剂 +5 February_Feb 2026-03-19 5/250 2026-03-19 23:51 by 23Postgrad
[考博] 申博26年 +3 八6八68 2026-03-19 3/150 2026-03-19 19:43 by nxgogo
[考研] 一志愿西安交通大学材料工程专业 282分求调剂 +5 枫桥ZL 2026-03-18 7/350 2026-03-19 14:52 by 功夫疯狂
[考研] 考研求调剂 +3 橘颂. 2026-03-17 4/200 2026-03-17 21:43 by 有只狸奴
[考研] 277调剂 +5 自由煎饼果子 2026-03-16 6/300 2026-03-17 19:26 by 李leezz
[考研] 308求调剂 +4 是Lupa啊 2026-03-16 4/200 2026-03-17 17:12 by ruiyingmiao
[考研] 302求调剂 +4 小贾同学123 2026-03-15 8/400 2026-03-17 10:33 by 小贾同学123
[考研] [导师推荐]西南科技大学国防/材料导师推荐 +3 尖角小荷 2026-03-16 6/300 2026-03-16 23:21 by 尖角小荷
[考研] 304求调剂 +4 ahbd 2026-03-14 4/200 2026-03-16 16:48 by 我的船我的海
[考研] 327求调剂 +6 拾光任染 2026-03-15 11/550 2026-03-15 22:47 by 拾光任染
信息提示
请填处理意见