版塊導(dǎo)航: 正在加載中...

應(yīng)《網(wǎng)絡(luò)安全法》要求，自2017年10月1日起，未進(jìn)行實(shí)名認(rèn)證將不得使用互聯(lián)網(wǎng)跟帖服務(wù)。為保障您的帳號(hào)能夠正常使用，請(qǐng)盡快對(duì)帳號(hào)進(jìn)行手機(jī)號(hào)驗(yàn)證，感謝您的理解與支持！

24小時(shí)熱門版塊排行榜

返回列表

glazio

鐵蟲 (小有名氣)

應(yīng)助: 1 (幼兒園)
金幣: 1023.4
散金: 12
帖子: 124
在線: 74.8小時(shí)
蟲號(hào): 453798
注冊(cè): 2007-11-07
專業(yè): 材料

[求助] 用python求7個(gè)原子的回轉(zhuǎn)半徑

寫代碼上我是門外漢，我自己算出來的結(jié)果懷疑不對(duì)，在這里請(qǐng)python達(dá)人指點(diǎn)一下，能給出你們的計(jì)算腳本讓我拜讀是最好了，呵呵。

1.  回轉(zhuǎn)半徑網(wǎng)頁請(qǐng)見這里

2.  代碼要求：因?yàn)檫@個(gè)腳本是在一個(gè)模擬軟件中使用的，該軟件使用基本的python環(huán)境而沒有scipy、numpy等模塊，所以在腳本里只能用python內(nèi)置的一些函數(shù)和命令

3.  七個(gè)原子的三維坐標(biāo)和我做的代碼如下

def f_Rg(list):
n=len(list)
sum_x=sum_y=sum_z=0
for list_sub1 in list:
      sum_x+=list_sub1[0]
      sum_y+=list_sub1[1]
      sum_z+=list_sub1[2]
rc=[1.0*sum_x/n,1.0*sum_y/n,1.0*sum_z/n]  #center of gravity, assuming all atoms have equal weight

sum_Rg2=0
for list_sub2 in list:
      sum_Rg2+=(list_sub2[0]-rc[0])**2+(list_sub2[1]-rc[1])**2+(list_sub2[2]-rc[2])**2

result=[rc,(1.0*sum_Rg2/n)**0.5]
return result

mol = [ [9.6397, 10.8355, 5.5044], \
         [7.6624, 11.6655, 5.5672], \
         [9.1239, 10.4822, 4.5657], \
         [8.1561, 10.4279, 4.4660], \
         [7.3390, 10.9163, 4.8606], \
         [7.4566, 12.0949, 4.5525], \
         [7.9986, 12.6564, 5.3909]]

Gy = f_Rg(mol)
print format( 'Center=(%.2f, %.2f, %.2f),  Rg=%.2f' %(Gy[0][0], Gy[0][1], Gy[0][2], Gy[1]) )

4.  以上腳本的計(jì)算結(jié)果如下
Center=(8.20, 11.30, 4.99),  Rg=1.22
下圖中綠色球就是7個(gè)原子，紅色箭頭起始點(diǎn)是重心rc，箭頭長(zhǎng)度就是Rg，怎么看怎么覺得計(jì)算結(jié)果是錯(cuò)的

7 atoms

[ Last edited by glazio on 2012-6-21 at 13:00 ]

回復(fù)此樓

» 猜你喜歡

279分求調(diào)劑一志愿211 已經(jīng)有16人回復(fù)
求調(diào)劑一志愿武漢理工大學(xué)材料工程（085601）已經(jīng)有5人回復(fù)
上海電力大學(xué)材料防護(hù)與新材料重點(diǎn)實(shí)驗(yàn)室招收調(diào)劑研究生（材料、化學(xué)、電化學(xué)，環(huán)境）已經(jīng)有3人回復(fù)
急發(fā)核心期刊論文已經(jīng)有5人回復(fù)
291 求調(diào)劑已經(jīng)有5人回復(fù)
298-一志愿中國(guó)農(nóng)業(yè)大學(xué)-求調(diào)劑已經(jīng)有11人回復(fù)
一志愿南京理工大學(xué)085701資源與環(huán)境302分求調(diào)劑已經(jīng)有7人回復(fù)
316求調(diào)劑已經(jīng)有7人回復(fù)
材料與化工考研調(diào)劑已經(jīng)有4人回復(fù)
一志愿重慶大學(xué)085700資源與環(huán)境，總分308求調(diào)劑已經(jīng)有7人回復(fù)

» 本主題相關(guān)價(jià)值貼推薦，對(duì)您同樣有幫助:

通過蛋白質(zhì)骨架還原全原子坐標(biāo) 已經(jīng)有6人回復(fù)
兩原子間出現(xiàn)(3,-3)關(guān)鍵點(diǎn)和兩個(gè)（3，-1）關(guān)鍵點(diǎn)，如何處理？已經(jīng)有4人回復(fù)
誰有各原子軌道半徑數(shù)據(jù)+如何確定某元素的自旋多重度？已經(jīng)有3人回復(fù)
原子半徑隨壓強(qiáng)增大會(huì)減小嗎已經(jīng)有4人回復(fù)
求元素的原子實(shí)半徑表，或者原子實(shí)半徑的求法已經(jīng)有5人回復(fù)
RAW264.7細(xì)胞長(zhǎng)得很慢嗎？已經(jīng)有12人回復(fù)
原子間距小于贗勢(shì)的截?cái)喟霃綍r(shí)的情況應(yīng)如何處理呢？已經(jīng)有6人回復(fù)
【求助】現(xiàn)在大半徑碳納米管的摻雜建立超胞要計(jì)算多少個(gè)原子已經(jīng)有5人回復(fù)
【求助】請(qǐng)問 xrd 和離子半徑有關(guān)，還是原子核半徑有關(guān)？已經(jīng)有11人回復(fù)
【求助】Python 中求解微分方程 odeint 的用法已經(jīng)有11人回復(fù)
【求助】Abinit中計(jì)算PDOS時(shí)怎樣確定投影原子半徑？已經(jīng)有4人回復(fù)
【請(qǐng)教】diamond填充圖的原子半徑怎么調(diào)小一些已經(jīng)有4人回復(fù)
【求助】氧化鉍熱膨脹系數(shù) 已經(jīng)有5人回復(fù)
【求助】請(qǐng)問從哪里可以查到元素的原子半徑？可以表明參考文獻(xiàn)的那種已經(jīng)有12人回復(fù)
【分享】世界最小晶體管問世：僅由7個(gè)原子構(gòu)成已經(jīng)有5人回復(fù)
【素材】原子(離子)半徑(pm)周期表已經(jīng)有131人回復(fù)

1樓 2012-06-21 12:31:12

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

libralibra

至尊木蟲 (著名寫手)

驃騎將軍

程序強(qiáng)帖: 40
應(yīng)助: 817 (博后)
金幣: 12914.1
紅花: 64
帖子: 2238
在線: 287.3小時(shí)
蟲號(hào): 696514
注冊(cè): 2009-02-05
專業(yè): 計(jì)算機(jī)軟件

如果按照這個(gè)公式來計(jì)算,的出來的是個(gè)很小的數(shù)
你的答案跟這個(gè)中始終的r_mean是一致的,也就是the mean position of the monomers
不過真心不懂不要做的東西,光能根據(jù)公式看看

CODE:

#! /usr/bin/env python

#

mol = [[9.6397, 10.8355, 5.5044], \

            [7.6624, 11.6655, 5.5672], \

            [9.1239, 10.4822, 4.5657], \

            [8.1561, 10.4279, 4.4660], \

            [7.3390, 10.9163, 4.8606], \

            [7.4566, 12.0949, 4.5525], \

            [7.9986, 12.6564, 5.3909]]

# r_mean

rmean = [x/len(mol) for x in reduce(lambda x,y:[x[i]+y[i] for i in range(len(x))],mol)]

# r_k - r_mean

subsum = [map(lambda x,y: x-y, mol_item,rmean) for mol_item in mol]

# sum(r_k - r_mean)^2

sumsqure = map(lambda x: [pow(y,2) for y in x], subsum)

# 1/N * sum(r_k - r_mean)^2

Rg_2 = [x/len(sumsqure) for x in reduce(lambda x,y:[x[i]+y[i] for i in range(len(x))],sumsqure)]

# sqrt(1/N * sum(r_k - r_mean)^2

Rg = [pow(x,0.5) for x in Rg_2]

# print result

print 'r_mean = ', rmean

print 'Rg^2 = ', Rg_2

print 'Rg = ', Rg

結(jié)果

CODE:

>>> ================================ RESTART ================================

>>> 

r_mean =  [8.196614285714286, 11.296957142857142, 4.986757142857144]

Rg^2 =  [0.650244935510204, 0.6282219424489798, 0.20302449959183683]

Rg =  [0.8063776630774219, 0.7926045309288736, 0.4505824004461746]

>>>

贊一下

回復(fù)此樓

matlab/VB/python/c++/Java寫程序請(qǐng)發(fā)QQ郵件:790404545@qq.com

2樓2012-06-21 16:41:30

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

libralibra

至尊木蟲 (著名寫手)

驃騎將軍

程序強(qiáng)帖: 40
應(yīng)助: 817 (博后)
金幣: 12914.1
紅花: 64
帖子: 2238
在線: 287.3小時(shí)
蟲號(hào): 696514
注冊(cè): 2009-02-05
專業(yè): 計(jì)算機(jī)軟件

【答案】應(yīng)助回帖

★ ★ ★ ★ ★
感謝參與，應(yīng)助指數(shù) +1
glazio: 金幣+5, ★★★很有幫助 2012-06-21 18:31:45

CODE:

#! /usr/bin/env python

#

mol = [[9.6397, 10.8355, 5.5044], \

            [7.6624, 11.6655, 5.5672], \

            [9.1239, 10.4822, 4.5657], \

            [8.1561, 10.4279, 4.4660], \

            [7.3390, 10.9163, 4.8606], \

            [7.4566, 12.0949, 4.5525], \

            [7.9986, 12.6564, 5.3909]]

# r_mean

rmean = [x/len(mol) for x in reduce(lambda x,y:[x[i]+y[i] for i in range(len(x))],mol)]

# r_k - r_mean

subsum = [map(lambda x,y: x-y, mol_item,rmean) for mol_item in mol]

# sum(r_k - r_mean)^2

sumsqure = map(lambda x: [pow(y,2) for y in x], subsum)

# 1/N * sum(r_k - r_mean)^2

Rg_2 = [x/len(sumsqure) for x in reduce(lambda x,y:[x[i]+y[i] for i in range(len(x))],sumsqure)]

# sqrt(1/N * sum(r_k - r_mean)^2

Rg = [pow(x,0.5) for x in Rg_2]

# print result

print 'r_mean = ', rmean

print 'Rg^2 = ', Rg_2

print 'Rg = ', Rg

結(jié)果

CODE:

>>> ================================ RESTART ================================

>>> 

r_mean =  [8.196614285714286, 11.296957142857142, 4.986757142857144]

Rg^2 =  [0.650244935510204, 0.6282219424489798, 0.20302449959183683]

Rg =  [0.8063776630774219, 0.7926045309288736, 0.4505824004461746]

>>>

贊一下

回復(fù)此樓

matlab/VB/python/c++/Java寫程序請(qǐng)發(fā)QQ郵件:790404545@qq.com

3樓2012-06-21 16:42:31

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

glazio

鐵蟲 (小有名氣)

應(yīng)助: 1 (幼兒園)
金幣: 1023.4
散金: 12
帖子: 124
在線: 74.8小時(shí)
蟲號(hào): 453798
注冊(cè): 2007-11-07
專業(yè): 材料

感謝libra回復(fù)！代碼正在拜讀，其中map、lambda的方法都是我還沒學(xué)到的東西，需要去翻教科書來理解。

疑問在于，回轉(zhuǎn)半徑應(yīng)該是標(biāo)量。
例如矢量r的平方，其值等于r的模( |r| )的平方，即 r^2 = |r|^2

請(qǐng)libra有時(shí)間再解答一下這個(gè)疑問，謝謝。

贊一下

回復(fù)此樓

4樓2012-06-21 17:12:31

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

libralibra

至尊木蟲 (著名寫手)

驃騎將軍

程序強(qiáng)帖: 40
應(yīng)助: 817 (博后)
金幣: 12914.1
紅花: 64
帖子: 2238
在線: 287.3小時(shí)
蟲號(hào): 696514
注冊(cè): 2009-02-05
專業(yè): 計(jì)算機(jī)軟件

引用回帖:

4樓: Originally posted by glazio at 2012-06-21 17:12:31
感謝libra回復(fù)！代碼正在拜讀，其中map、lambda的方法都是我還沒學(xué)到的東西，需要去翻教科書來理解。

疑問在于，回轉(zhuǎn)半徑應(yīng)該是標(biāo)量。
例如矢量r的平方，其值等于r的模( |r| )的平方，即 r^2 = |r|^2

請(qǐng)libr ...

那可能是我計(jì)算錯(cuò)誤,我計(jì)算的時(shí)候
r =[x,y,z],r_mean=[mean(x),mean(y),mean(z(],然后(r-r_mean)的平方我直接是[x^2,y^2,z^2]了

贊一下

回復(fù)此樓

matlab/VB/python/c++/Java寫程序請(qǐng)發(fā)QQ郵件:790404545@qq.com

5樓2012-06-21 18:39:29

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

libralibra

至尊木蟲 (著名寫手)

驃騎將軍

程序強(qiáng)帖: 40
應(yīng)助: 817 (博后)
金幣: 12914.1
紅花: 64
帖子: 2238
在線: 287.3小時(shí)
蟲號(hào): 696514
注冊(cè): 2009-02-05
專業(yè): 計(jì)算機(jī)軟件

內(nèi)容已刪除

回復(fù)此樓

matlab/VB/python/c++/Java寫程序請(qǐng)發(fā)QQ郵件:790404545@qq.com

6樓2012-06-21 18:46:22

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

glazio

鐵蟲 (小有名氣)

應(yīng)助: 1 (幼兒園)
金幣: 1023.4
散金: 12
帖子: 124
在線: 74.8小時(shí)
蟲號(hào): 453798
注冊(cè): 2007-11-07
專業(yè): 材料

五體投地矣，拜謝！

回復(fù)此樓

7樓2012-06-21 18:55:07

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 glazio 的主題更新

返回列表

最具人氣熱帖推薦 [查看全部]		作者	回/看	最后發(fā)表

[考研] 328求調(diào)劑 +4	LHHL66 2026-03-23	4/200	2026-03-23 14:55 by lbsjt
[考研] 一志愿上海交大生物與醫(yī)藥專碩324分，求調(diào)劑 +4	jiajunX 2026-03-22	4/200	2026-03-23 14:21 by xiaotangjing
[考研] 0703 調(diào)劑 +3	我可以上岸的對(duì)?/a> 2026-03-16	6/300	2026-03-23 08:52 by 醉在風(fēng)里
[考研] 293求調(diào)劑 +3	濤濤Wjt 2026-03-22	5/250	2026-03-22 22:21 by jiangpengfei
[考研] 298求調(diào)劑一志愿211 +3	上岸6666@ 2026-03-20	3/150	2026-03-22 15:50 by ColorlessPI
[考研] 384求調(diào)劑 +3	子系博 2026-03-22	4/200	2026-03-22 11:04 by 搏擊518
[考研] 286求調(diào)劑 +10	Faune 2026-03-21	10/500	2026-03-21 23:34 by 314126402
[考研] 280求調(diào)劑 +11	咕嚕曉曉 2026-03-18	12/600	2026-03-21 22:40 by ACS Nano——
[考研] 初試 317 +7	半拉月丙 2026-03-20	7/350	2026-03-21 22:26 by peike
[考研] 材料工程專碩 348分求調(diào)劑 +3	冬辭. 2026-03-17	5/250	2026-03-21 18:47 by 學(xué)員8dgXkO
[考研] 0805材料320求調(diào)劑 +3	深海物語 2026-03-20	3/150	2026-03-21 15:46 by 無際的草原
[考研] 南昌大學(xué)材料專碩311分求調(diào)劑 +6	77chaselx 2026-03-20	6/300	2026-03-21 07:24 by JourneyLucky
[考研] 301求調(diào)劑 +10	yy要上岸呀 2026-03-17	10/500	2026-03-21 03:14 by JourneyLucky
[考研] 332求調(diào)劑 +4	ydfyh 2026-03-17	4/200	2026-03-21 02:20 by JourneyLucky
[考研] 一志愿重慶大學(xué)085700資源與環(huán)境專碩，總分308求調(diào)劑 +3	墨墨漠 2026-03-18	3/150	2026-03-21 00:39 by JourneyLucky
[考研] 274求調(diào)劑 +10	S.H1 2026-03-18	10/500	2026-03-20 23:51 by JourneyLucky
[考研] 308求調(diào)劑 +3	阿姐阿姐家啊 2026-03-18	3/150	2026-03-20 23:24 by JourneyLucky
[考研] 290求調(diào)劑 +7	^O^乜 2026-03-19	7/350	2026-03-20 21:43 by JourneyLucky
[考研] 295復(fù)試調(diào)劑 +8	簡(jiǎn)木ChuFront 2026-03-19	8/400	2026-03-20 20:44 by zhukairuo
[考博] 申博26年 +3	八6八68 2026-03-19	3/150	2026-03-19 19:43 by nxgogo

24小時(shí)熱門版塊排行榜

[求助] 用python求7個(gè)原子的回轉(zhuǎn)半徑

» 猜你喜歡

» 本主題相關(guān)價(jià)值貼推薦，對(duì)您同樣有幫助:

【答案】應(yīng)助回帖

» 本主題相關(guān)價(jià)值貼推薦，對(duì)您同樣有幫助: