亭亭五月天在线观看,亭亭五月天在线观看,国产最新av一区二区,国产 高清 中文字幕,99re热久久亚洲综合精品成人,熟妇 一区二区三区,一级做a爰片性色毛片武则天,美女的骚穴视频播放,国产美女午夜免费视频

24小時熱門版塊排行榜    

查看: 16051  |  回復: 63
【有獎交流】積極回復本帖子,參與交流,就有機會分得作者 luyaobao 的 765 個金幣 ,回帖就立即獲得 5 個金幣,每人有 1 次機會

luyaobao

木蟲 (著名寫手)


[交流] 分子動力學模擬及其LAMMPS實現

(課程尚未完成,持續(xù)更新中)(三萬字超長文預警)

前言

1. 分子動力學基礎

1.1 原理.

1.2 勢能

1.3 溫度、壓強

2. linux系統常見命令和操作

3. lammps介紹和并行版安裝(linux版)

4. lammps模擬基本流程

5. 建模和后處理軟件概述

6. 出錯解決思路

7. 專題一:納米流動(部分內容)

8. 專題二:納米流體(待更新)

9. 專題三:潤濕(部分內容)

9.1建模

10. 專題四:拉伸、壓縮、彎曲、扭轉、剪切(待更新)

11. 專題五:傳熱(待更新)

12. 專題六:粘度和擴散(待更新)

13. 專題八:摩擦和磨損(待更新)

14. 專題九:反應分子動力學reaxff(待更新)

15. 專題十:耗散粒子動力學(待更新)

16. 專題十一:自由能(待更新)

17. 專題十二:原子沉積和鍍膜(待更新)

前言

分子動力學模擬在科研中的使用越來越廣泛。lammps作為分子動力學中一款優(yōu)秀的開源軟件,使用率非常高。大量高水平的論文都是采用lammps完成的,其中不乏nature,science,prl,jacs,nature communication,pnas等頂刊。但是對于零基礎的同學,特別是本科生和低年級碩士生,lammps的學習曲線還是比較陡峭的。很多同學都是課題組內第一個使用分子動力學模擬和lammps,其困難可想而知。如果能有一個針對lammps的系列教程對學生掌握這個工具會大有裨益。本教程就是希望能夠在學生使用lammps初期降低學習難度,少走彎路,快速上手,將自己的精力集中在具體的科學問題上,而不是軟件的學習上。本教程假設你之前對linux,分子動力學和lammps一無所知,從頭開始講,因此比較淺顯,并未就某些問題進行深入介紹(其實是深入的我也不會,哈哈)。

我將分子動力學研究科學問題分為三個部分:軟件使用,分子動力學理論,科學問題發(fā)現。如果把分子動力學比作武林,那么軟件的使用僅僅是招式,分子動力學理論則是內功心法,而科學問題則是對手。要想成為一名武林高手,在江湖上闖出一番事業(yè),既要熟練掌握各種招式,也要有很深的內功修為,同時也要有一些值得尊敬的對手。只有用深厚的內功催動精妙的招式,才能打敗可敬的對手,受到江湖傳頌。本教程只是lammps的使用教程,它會教給你使用lammps的一招一式,但要將招式練得熟練,還需要自己在科研中積累大量模擬經驗。而內功心法的修煉則需要長時間的學習經典教材和文獻。好的對手必須要廣泛閱讀文獻,與導師深入溝通確定具有很好學術價值的科學問題。那么有沒有一種方法像武俠小說中的主角,掉下懸崖然后打通任督二脈,迅速成為高手。答案是當然是沒有的。但是作為物理學分支,我推薦閱讀《費曼物理講義》第一卷,可以起到類似的作用。總之,科研沒有坦途,靜心學習,努力思考必有所成。

“l(fā)ammps即large-scale atomic/molecular massively parallel simulator,可以翻譯為大規(guī)模原子分子并行模擬器,主要用于分子動力學相關的一些計算和模擬工作,一般來講,分子動力學所涉及到的領域,lammps代碼也都涉及到了。lammps由美國sandia國家實驗室開發(fā),以gpl license發(fā)布,即開放源代碼且可以免費獲取使用,這意味著使用者可以根據自己的需要自行修改源代碼。lammps可以支持包括氣態(tài),液態(tài)或者固態(tài)相形態(tài)下、各種系綜下、百萬級的原子分子體系,并提供支持多種勢函數。且lammps有良好的并行擴展性!薄园俣劝倏。lammps的官網是https://www.lammps.org/。在lammps的官網中有大量關于軟件的信息,可以好好探索一番。有幾個地方可以經?纯础5谝皇莗ublication部分https://www.lammps.org/papers.html,記錄引用lammps的所有文章,多看看對自己的研究有所啟發(fā)。maillist部分https://www.lammps.org/mail.html是lammps的支持社區(qū)。當你有問題的時候可以檢索maillist,你遇到的問題別人也遇到過。所以查看開發(fā)者對這些問題的解答。

如何學習lammps?你可以把本教程當做一個開始。lammps的手冊中詳細介紹了軟件的各個方面。一定要好好閱讀手冊。手冊的前四章要認真閱讀。經常使用的命令也要仔細閱讀。當你要實現某種功能的時候就把手冊打開看看命令。仔細閱讀學習lammps官方手冊是成為lammps高手的必經之路,閱讀學習一本分子動力學模擬經典教材是增加修為的關鍵。希望本教程能教會你lammps的基本招式。祝你好運!

lammps官網推薦了幾本分子動力學模擬教材:

books about molecular dynamics generally or lammps specifically

note that these are not endorsements of particular books. we simply want to make the lammps user community aware of them as potentially useful resources.

general md books:

allen & tildesley - computer simulation of liquids
frenkel & smit - understanding molecular simulation: from algorithms to application
griebel, knapek, zumbusch - numerical simulation in molecular dynamics: numerics, algorithms, parallelization, applications
tuckerman - statistical mechanics: theory and molecular simulation
books about lammps specifically:

mubin & li - extending and modifying lammps (see below)

1. 分子動力學基礎

1.1 原理

本節(jié)只是簡要介紹分子動力學的原理,如果要深入學習相關內容,可閱讀前言中推薦的教材。

分子動力學的基礎是牛頓力學,也即經典力學。經典力學中有三個主要內容:質點、力和運動。牛頓第二定律是經典力學的核心,下面兩個質點的方程大家一定不陌生。


我們把他寫成離散的形式,也稱為差分形式:




如果delta_t足夠小那么離散形式可以以足夠的精度近似連續(xù)形式。首先我們假設質點的受力之依賴于它的位置,也就是說只是位置的函數。事實上這是一個很普遍成立的假設,比如重力或者萬有引力,在質量確定的情況下只依賴于兩個物體之間的相對距離。在比如靜電力,在電荷確定的情況下只依賴于兩個帶電質點的相對距離。有了這個假設,我們考慮一個一維的簡單情形。有兩個質點放置在光滑的地面上,如下圖。


在t=0時刻,質量為m_a位于x_0^a的質點a以速度v_0^a朝質量為m_b位于x_0^b且靜止(即v_0^b=0)的質點b運動。兩個質點之間的受力f_ab只是質點a與質點b之間距離x_ab的函數,也即只要知道x_ab,我們就能計算出來f_ab。然后,我們要問的問題是在后續(xù)的時間中兩個質點是怎么運動的?我們這時設定一個很小的delta_t。我們知道初始速度那么有



這樣我們就知道了t_0+delta_t時刻,質點a和b的位置和速度。由新的位置又可以得到新的受力。重復上述公式就可以知道后續(xù)所有時刻兩個質點的位置和速度,也就是知道兩個質點的運動細節(jié)?偨Y一下我們結算所有時刻質點的運動細節(jié)要提前知道哪些信息。質點的質量、初始速度、初始位置和依賴于位置的受力函數。如果我們再選一個足夠小的delta_t,那么我們就能夠準確的得到質點a和b的所有運動信息。

我們把情形擴展到二維。在一個光滑的臺球桌面上,白球在初始時刻以確定的速度撞向其他擺好的臺球。如果我們把臺球都看成質點,又知道所有臺球的質量、初始速度、初始位置和依賴于位置的受力函數,我們就能預測擊球后所有時刻臺球的運動軌跡,從而判斷臺球是否能夠進袋。我們再把情形擴展到三維?紤]宇宙中只存在太陽和太陽系的八大行星。我們是否能夠通過萬有引力去預測所有行星的軌跡。答案當然是可以的。只是此時我們知道某個行星除了來自太陽的引力,還有來自其他行星的應力。我們在計算某個行星的受力時,要考慮該行星與其余所有天體之間的受力。這就是經典力學能夠告訴我們的東西。

在這個概念上,拉普拉斯提出了一個想法。如果存在一只神奇的動物,它能基于經典力學解算出世界上所有東西的位置和速度并能加以分析,那么世界上發(fā)生的一切及其背后的原因都一覽無余的呈現在它眼前。這就是經典力學的世界觀和認識論。分子動力學的理論就是基于這樣的認識方法。像上面的例子一樣,分子動力學的目的就是為了求解任意時刻所有組成物體原子的位置和速度,從而認識物質世界。


為什么這個方法具有很大的威力呢?這是多年科學實踐告訴我們的。事實上著名物理學家理查德·費曼對這個問題有著精彩的描述。理查德·費曼在其傳世物理教材《費曼物理講義》中提出這樣一個有意思的問題:“假如由于某種大災難,所有的科學知識都丟失了,只有一句話傳給下一代,那么怎樣才能用最少的詞匯來表達最多的信息呢?”費曼認為這句話應該是原子假設:“所有物體都是用原子構成的——這些原子是一些小小的粒子,它們一直不停地運動著,當彼此略微離開地時互相吸引,當彼此過于靠近時又相互排斥!边@句話之所以重要是因為基于這個原子假設我們能解釋巨量自然界中的現象。比如,氣體壓強來自于大量氣體分子對壁面撞擊;水蒸發(fā)吸熱來自于水分子從體相水中逃逸至外界環(huán)境而帶走能量所致。。。。。更多詳細例子可閱讀《費曼物理講義》第一卷的相關內容。

費曼把原子假設放在如此高度的原因總結起來就是——掌握組成物體內部原子的行為就可解釋和理解物體的性質和行為。只是這個時候我們關注的是大量原子在一起表現出來的整體特性,而不是關注具體某個原子的運動細節(jié)。事實上,實踐告訴我們大量原子在一起就會表現出某些特定的確定的性質和行為,也就我們實際生活生產中能直接測量和感受的性質和行為。到這里我們可以看到分子動力學實際上就是一只拉普拉斯獸。分子動力學的基本任務就是獲取物體在任意時刻組成原子的所有位置和動量然后利用統計力學知識理解物體的性質和行為。如果這樣還不能說服你分子動力學的實力。就講一個事實——2013年諾貝爾化學獎頒給了創(chuàng)立跨尺度分子模擬的三位科學家。通過分子動力學模擬我們可以深刻理解我們所研究的對象。目前各種硬件超強的計算能力為分子動力學創(chuàng)建了一個大有可為的時代!!

分子動力學的基本任務就是獲得研究對象不同時刻的位置和動量,然后基于統計力學知識獲得想要的物理量,解釋對象的性質和行為。因此,分子動力學的模擬流程超級簡單。第一步,設置研究對象組成原子或者粒子(下面統一稱為粒子)初始位置和速度;第二步,基于粒子的位置計算每個粒子受到的合力,并基于牛二定律計算粒子的加速度;第三步,基于加速度計算粒子下一時刻的速度;第四步,基于下一時刻的速度計算下一時刻的位置;第五步,循環(huán)第二步到第四步的過程,得到一系列時刻粒子的位置和速度;第六步,基于位置和速度信息得到描述對象性質和行為的物理量。


分子動力學的基本原理就是這么簡單,正因為簡單所以有效。需要指出的是在實際編程求解的過程考慮到計算精度等原因上述過程在某些細節(jié)上會略有改動。但是,基本上你可以認為這就是分子動力學的原理和流程。上述過程由lammps幫你實現,你需要做的就是設置研究對象的基本信息,告訴給lammps,然后開始模擬。

1.2 勢能

前面我們知道要獲得所有原子的軌跡,需要知道原子的質量,初始位置,初始速度和依賴于位置的受力函數。原子的質量,初始速度和初始位置都是我們人為設定的值。不同研究對象的原子質量,初始速度和位置很容易獲得。不同研究對象更本質的差別在于依賴于位置的受力函數,這個函數又被稱為相互作用。相互作用的形式和參數決定了研究對象的性質和行為。在分子動力學中,這個相互作用通過原子之間的勢能進行描述,因此又稱為勢能函數或勢函數。而受力是勢能對位置的負梯度,因為實踐表明原子之間的受力都是保守力。舉個例子,重力和重力勢能。取豎直向上為正方向,則重力表達式為


重力勢能的表達式是


即重力勢能對位置的負梯度就是重力


為什么要用勢能,而不是受力來表達相互作用呢?這是因為在絕大多數模擬中采用的是笛卡爾坐標。勢能是標量,而力是矢量,數學上處理起來更容易。在計算受力的時候,直接用勢能對兩個原子之間的相對位置進行求導,然后將這個力乘以方向余弦就能得到三個方向上的受力。舉例說明,在三維空中有兩個原子


兩個原子之間的距離


那么原子對原子的受力就為


這樣不管是程序的簡潔性還是計算效率都會很好。所以要開展分子動力學計算,首先必須確定所有原子之間的勢能函數,即勢能關于兩個原子之間相對位置的函數。確定函數分為兩個部分:函數表達式和表達式中的常數值。對于一個對象,描述其相互作用的勢能函數表達式和表達式中常數值的集合稱為該對象的力場。通常一種力場可以用來描述一類對象。分子動力學模擬一般研究的對象可以分為:惰性氣體,分子物質,金屬晶體和非金屬晶體。分子動力學模擬絕大多數情況用來研究固體和液體,即凝聚態(tài)。對于氣體分子動力學模擬效率太低,有更好的辦法進行研究。簡單液體原子之間相互作用只具有范德瓦爾斯作用的物質。描述范德瓦爾斯作用最常用的表達式是lennard-jones12/6勢能。通常用lj12/6描述的物質也稱為簡單液體。





上圖是lj/126勢函數的圖像。對于該勢函數有幾個重要的點要說一下。勢能的最低點的到零點之間的絕對值為ε,稱為勢能函數的特征能量,勢能的零點位于處r=σ,稱為勢能函數的特征長度。勢能函數的最低點位于處r=2^1/6σ,在該點處勢能對位置的導數是零,所以該點處受力也為零。該點也稱為平衡位置,我們可以看到在平衡位置出勢能是負的。所以當體系處于平衡狀態(tài)時,原子之間的相對距離大概也在平衡位置附近,所以勢能為負的。這也就是為什么大多數情況分子動力學算出來的勢能都是負的。在r<2^1/6σ,兩個原子之間表現為相互排斥,在r>2^1/6σ,兩個原子之間表現為相互吸引。從勢函數的圖像上可以看出在較大的時候,勢能趨于零,其導數也趨于零。也就是當兩個原子距離很遠的時候,原子之間的相互作用可以忽略不計。因此在實際計算的時候都會使用截斷的lj/126。截斷距離成為截斷半徑,其值大概在。就像重力勢能一樣,勢能的計算需要一個參考位置。lj/126的參考位置選擇在了無限遠處。類似lj/126勢能以及其他拓展的lj勢能,都稱為對勢,其含義是兩個原子之間的相互作用只由這一對原子決定,與其它相鄰原子沒有關系。分子物質原子之間有化學鍵連接的物質。比如常見的水,甲烷,乙醇,醋酸,氯化鈉,聚乙烯,酯類,氨基酸,rna,dna等。分子物質是分子動力學模擬應用最廣的對象之一。為了描述分子物質研究者開發(fā)出了一系列力場。值得說明的是力場是一個經驗性的東西。研究者通過自己的物理直覺和反復試驗,最終建立起適合某一類物質的力場。對于分子物質,相互作用包括兩類:非鍵結勢能和鍵結勢能。非鍵結用來描述沒有化學鍵連接原子之間的相互作用,鍵結勢能用來描述有化學鍵連接原子之間的相互作用,如下圖。


在一個體系中任意兩個非鍵結的原子,如果其距離小于截斷半徑那么都要考慮它們之間的非鍵結相互作用。對于大多數力場非鍵結勢能用lj12/6描述。對于鍵結原子之間的相互作用則通過鍵結勢函數進行描述。鍵結相互作用分為鍵能(bond),鍵角能(angle),二面角能(dihedral),離平面能(improper)。兩個以化學鍵相連的兩個原子都具有鍵能相互作用,當化學鍵偏離平衡長度的時候兩個原子就會出現相互作用。三個化學鍵連接的相鄰原子之間均具有鍵角能,當三個原子之間的角度偏離平衡角度的時候就會產生鍵角能。四個化學鍵連接的相鄰原子之間會存在二面角能相互作用,當四個原子構成的二面角偏離平衡二面角的時候四個原子之間會出現二面角相互作用。四個相鄰化學鍵連接的相鄰原子之間有時還會存在離平面勢能作用,當第四個原子偏離前三個原子所構成平面的平均夾角時,四個原子之間會存在離平面能相互作用。不同的分子會具有不同的鍵結勢能。如氧氣分子就只具有鍵能相互作用,水分子具有鍵能和鍵角相互作用,十二烷分子具有鍵能,鍵角能和二面角相互作用分子。元素構成較為復雜的分子,如rna,就會同時具有四種鍵結相互作用。具體一個分子具有哪些鍵結能相互作用,由描述該分子的力場所定義。你只需要使用合適的軟件在建模的時候會自動定義好各種鍵結信息和勢能函數的參數。值得指出的是,在分子體系中原子一般都帶有電荷,所以非鍵結原子之間會存在靜電相互作用。在大多數力場中,原子的電荷都是部分電荷,也就是電荷一般都是零點幾幾幾。這是由于在共價鍵中任何一個原子都并沒有完全得到或失去一個電子,只是部分的得到或失去電子,因此在描述其電荷的時候都是零點幾幾幾。少數強電解質會有整數電荷,比如氯化鈉中的氯離子和鈉離子,分別是-1和+1電荷。這里的電荷單位是一個電子或質子所帶的電荷量,也就是單位電荷量。

為什么會采用多種鍵結勢能函數來描述分子體系呢?這是因為鍵結的原子之間電子云的相互作用會對分子的構型進行限制,導致分子會具有特定的幾何形狀。分子動力學為描述這種情形就采用這樣鍵能組合的形式。這是一種經驗的方法。實踐表明這是一種很好的辦法。只要有合適的勢函數形式和參數,采用的力場就能夠準確描述分子的構型和之間的相互作用。目前研究者針對分子體系已經開發(fā)出多套力場。
金屬晶體和非金屬晶體

這兩類物質具有特殊的結構,研究者為它們開發(fā)了額外的勢函數,稱為多體勢函數,來描述原子之間的作用。所謂多體就是說兩個原子之間的相互作用不僅取決于這兩個原子還與其周圍的原子有關系。針對某種物質研究者已經開發(fā)出相應的力場參數,只需要搜索一般就能找到對應的勢函數參數或文件。我對這些勢函數不太了解,就不在這里贅述了。

以上就是關于勢能和勢函數的內容。需要指出的是為了提高模擬精度,研究者還針對開發(fā)出了針對某一種物質的勢函數。當你要模擬這些物質時,直接查找原始文件中的參數就行。關于力場的選擇,當然是選針對要研究體系的力場。但是對于lammps來說建模是個很麻煩的事。如果你能找到一種力場描述的體系,并且成功建模,那你就偷著樂吧。不同力場之間的差別并不是那么明顯。事實上,分子動力學模擬最大的價值就在于給出研究對象的一般性規(guī)律和背后的機理,而不是計算出某些精確的數值。而大多數力場能夠捕捉研究對象的本質特征,這就基本上能夠保證這個力場能夠再現研究對象的本質規(guī)律,這對于我們來說就已經足夠了。當然,力場的開發(fā)一直在發(fā)展中以越來越精確的描述研究對象。比如幾年來發(fā)展出來的機器學習力場就能夠以量子的精度開展分子動力學模擬。這是分子動力學模擬方法中的一大進展。當你要驗證一個力場的時候,計算研究對象你所關心的物理量的數值和實驗或者文獻對比,如果接近,力場就能用。

1.3 溫度、壓強(這里寫的有點簡單了,先就這樣吧) 溫度對于分子動力學模擬是一個很重要的物理量,因為分子動力學一般用來模擬凝聚態(tài),如果溫度過高就氣化了。初中物理告訴我們溫度是原子熱運動劇烈程度的度量。原子熱運動越劇烈,溫度越高。也就是說溫度和原子的熱運動速度是直接相關的。根據統計力學三維系統中溫度和原子速度具有以下對應關系


其中,是n原子總個數,k_b是玻爾茲曼常數。通過調整原子的速度就可以實現控制系統的溫度。

壓強是另一個重要的物理量。壓強在分子動力學中這樣定義的


其中,v是系統體積。第一項來自原子熱運動的貢獻。第二項來自原子之間相互作用的貢獻,稱為維里項。通過調整系統的體積和原子的位置就可以實現控制系統的壓強。對于液體來說壓強就等于我們平時所接觸到的壓強。而對于固體來說這里的壓強實際上就是內應力。事實上,壓強其實就是系統內部應力張量的對角線分量。壓強和溫度都是強度量,就是說其值與系統規(guī)模無關。與之相對應的是廣延量,比如能量,當系統規(guī)模增大時,系統能量增加。

更多分子動力學的細節(jié)將結合lammps使用過程進行介紹。

2. linux系統常見命令和操作

lammps是在針對linux系統開發(fā)的。雖然windows系統也能用但是會出現很多問題,所以強烈推薦你使用linux系統運行l(wèi)ammps。要搞模擬你一定要熟悉linux的命令行操作邏輯,要習慣使用腳本來寫lammps命令。很多軟件包括開源的和商業(yè)的都是采用腳本控制的方式來執(zhí)行的。

在linux系統中使用最多的不是鼠標點點點,而是在終端中輸入命令行來執(zhí)行操作。執(zhí)行一條命令是指在終端中輸入該命令,然后回車執(zhí)行。在ubuntu系統中可以使用ctrl+alt+t來打開終端,其他linux操作系統可以在軟件列表中點擊打開終端。這里只列舉最常用的命令。

cd命令用來進入某個目錄,比如cd /home就是進入home文件夾。在linux中.表示當前目錄,..表示上一級目錄。ls命令列出當前目錄下的所有文件和文件夾。pwd命令返回當前目錄的路徑。mkdir創(chuàng)建文件夾,如mkdir software表示在當前目錄下創(chuàng)建software文件夾,mkdir /home/softeare,表示在/home文件夾下創(chuàng)建software文件夾。rm命令刪除文件或文件夾。tar解壓縮命令當壓縮文件的擴展名是.tar.gz時可以使用tar -xvzf 壓縮文件名,將壓縮文件解壓。which命令,查找可執(zhí)行文件的路徑,如which lmp_g++_mpich就是查找可執(zhí)行文件lmp_g++_mpich的所在路徑,如果返回路徑那么這個命令就可以被系統找到,如果什么也沒返回就是找不到。echo命令返回變量的值,如echo $path就是返回path的值。環(huán)境變量是指系統或軟件定義的變量,一般用于文件查找。最常見的環(huán)境變量就是path,該變量記錄了眾多可執(zhí)行文件的位置。如果一個可執(zhí)行文件的路徑被添加到path中了,那么在運行該文件時可以直接在終端中寫該文件的名字,如果沒有被添加那么就要寫出該文件的絕對路徑才能執(zhí)行該文件。將一個路徑添加到path中需要修改.bashrc文件。輸入執(zhí)行vim ~/.bashrc,就可以打開.bashrc文件,按i鍵進入編輯模式,將光標定位在文件末尾并添加以下一行

path=$path:<路徑名>,按esc鍵退出編輯模式,輸入:wq后回車,這樣就保存好了剛才的編輯,并退出了.bashrc文件回到了終端。在終端輸入source ~/.bashrc就可以是剛才的操作發(fā)生作用。sudo是臨時取得root權限命令,當你要在系統目錄中添加或修改文件的時候就需要把sudo加到執(zhí)行命令前面,如sudo cp ./lmp /usr/bin就是把當前目錄下的lmp文件拷貝到系統目錄/usr/bin中。帶有sudo的命令回車執(zhí)行后需要輸入用戶密碼,當你在輸入密碼的時候終端上不會顯示密碼,但是確實是輸進去了,輸完回車如果密碼正確就開始執(zhí)行命令了。其他命令可以百度能夠找到大量資料。當你在使用linux系統時需要問一個問題,你就百度,一定能找到答案。

3. lammps介紹和并行版安裝(linux版)

這里介紹在ubuntu在編譯并行版lammps并執(zhí)行。首先假設你已經有了一臺安裝了ubuntu系統的電腦。你的用戶名叫me。電腦聯網。分別在lammps官網,mpich2官網和fftw官網下載對應的lammps,mpich和fftw的最新版源代碼,然后

第一步按ctrl+alt+t打開一個終端,輸入執(zhí)行以下命令

sudo apt-get update

sudo apt-get install build-essential

sudo apt-get install vim

第二步輸入執(zhí)行

mkdir /home/me/software

mkdir /home/me/software/mpich

mkdir /home/me/software/fftw

第三步將下載好的三個源代碼的壓縮文件拷貝進/home/me/software目錄下

第四步輸入執(zhí)行

tar -xvf <lammps壓縮文件的名稱>

tar -xvf <mpich壓縮文件的名稱>

tar -xvf <fftw壓縮文件的名稱>

第五步進入解壓出來的mpich的目錄執(zhí)行

./configure -prefix=/home/me/software/mpich --disable-fortran

make

make install

第六步進入解壓出來的fftw的目錄執(zhí)行

./configure -prefix=/home/me/software/fftw

make

make install

第七步輸入執(zhí)行

vim ~/.bashrc

打開.bashrc文件后按i鍵進入邊界模式,將光標定位在文件末尾并添加以下一行

path=$path:/home/me/software/mpich/bin

按esc鍵退出編輯模式,輸入:wq后回車,這樣就保存好了剛才的編輯,并退出了.bashrc文件回到了終端。在終端輸入

source ~/.bashrc

which mpirun

若輸出了mpirun的目錄則mpich安裝配置成功

第八步終端輸入執(zhí)行以下命令

cd /home/me/software/解壓出來lammps文件夾名稱/src

第九步打開目錄/home/me/software/解壓出來lammps文件夾名稱/src/make/options下的makefile.g++_mpich文件。在文件以下位置修改為

mpi_inc = -dmpich_skip_mpicxx -dompi_skip_mpicxx=1 -i/home/me/software/mpich/include

mpi_path = -l/home/me/software/mpich/lib

mpi_lib =

fft_inc = -i/home/me/software/fftw/include

fft_path = -l/home/me/software/fftw/lib

fft_lib =

保存退出。

第十步在終端依次執(zhí)行

make yes-molecule

make yes-kspace

make yes-rigid

make yes-manybody

make g++_mpich

等待編譯完成src目錄下會生成lmp_g++_mpich的可執(zhí)行文件,將此可執(zhí)行文件拷貝至你的in文件所在的位置輸入

mpirun -np <并行核數> ./lmp_g++_mpich -in <in文件名字>

運行正常則編譯成功。

4. lammps模擬基本流程

lammps的模擬需要寫一個腳本,稱之為in文件。我將一個典型的in文件書寫分為六個部分:1、定義模擬的基本規(guī)則,2、定義對象的幾何模型,3、定義對象的力場參數,4、定義分組信息,5、弛豫部分,6、正式模擬過程。

第一部分 定義模擬的基本規(guī)則

根據前面說的分子動力學原理,那么第一步就是設置粒子的初始位置和速度?等等。在這之前我們還需要做一些前期準備工作。首先是設置單位制,即涉及到的物理量的單位都是什么。為什么要這么做呢?就是為了能夠定量描述我們的對象;叵朐诟咧械闹麑嶒灐∏蚱綊。我們怎么描述小球的運動呢?我們會說小球經過多少秒,向前運動了多少米,下落了多少米。這里就采用了國際制單位描述小球的運動。只有先確定了單位制,我們才好定量描述研究對象。不同的研究對象,為了描述方便會采用不同的單位制。比如描述高鐵的速度我們會用km/h,而不是m/s。同樣,在lammps中面對不同的研究對象,lammps內置了不同的單位制。所以in文件中的

第一條命令要寫的就是設置單位制,即寫units命令。語法就是、

units <style> #舉例:units real

lammps所有的單位制可在手冊中查到。常見的單位系統有,描述分子體系的real,描述金屬體系的metal,描述簡單液體的lj。

第二條命令是設置邊模擬維度,即你模擬的是三維問題還是二維問題。通常我們的模擬都是三維的。語法是

dimension <n> #舉例:dimension 3

第三條命令是設置邊界條件,即boundary命令。通常我們希望描述對象的體相性質。但是,看看阿佛加德羅常數的值我們就可以知道實際的對象都會包含超級多的粒子。這么多的粒子的愛計算機無法承受。那怎么辦呢?分子動力學中采用了有限代替無限的做法。就是取對象的一部分然后加上周期性邊界條件來代替幾乎無限大的對象。如圖中所示所謂周期性邊界條件就是粒子如果從邊界穿出去,就會從對面穿回來。這樣就實現了以有限的粒子數代替了無限大的對象。由邊界所圍起來的區(qū)域通常稱為模擬盒子。模擬盒子的形狀最常用的就是長方體。之所以采用長方體是因為編程容易處理。


設置邊界條件在分子動力學模擬中優(yōu)先級很高,所以第二條命令要寫邊界條件,語法是

boundary <x> <y> <z> #舉例:boundary p p p

boundary p p p就是指三個方向上都是周期性邊界條件。這個設置在lammps中最常用。lammps中還提供了其他類型的邊界條件,常用的有固定邊界條件(f)和縮放邊界條件(s)。當我們需要在某個方向上不想采用周期性時就可用f或s。f作用就是粒子從邊界跑出去就回不來了,計算就把它丟掉了。s的作用是邊界的位置會隨著粒子的位置而縮放,當最靠近邊界的粒子往外移動時邊界位置也跟著移動,反之亦然。需要指出的是每個方向上都有兩個邊界,當寫一個參數是如p就表示該方向上兩個邊界都是周期性,f和s也如此。但是,還可以這么寫boundary p fs p,這就表示在x和z方向上是周期性邊界,在y方向上下面是f,上面是s。因此寫該條命令時根據需要進行設置。

第四條命令是設置建立鄰居列表的,即neighbor命令。根據分子動力學原理我們需要計算每個原子的受力。而任何一個粒子的受力都來自于其周圍粒子給它的合力。但是,我們知道只有距離較近的粒子才會對該粒子產生受力,距離較遠的粒子給與的力可以忽略不計。這就是鄰居的含義。所以要計算一個粒子所受的合力,首先要找到該粒子周圍有哪些鄰居,即建立鄰居列表。每一個粒子都有自己的鄰居列表。neighbor的語法是

neighbor <skin> <style> #舉例:neighbor 3 bin

所謂鄰居就是距離中心粒子一定范圍內的粒子,這個范圍稱之為作用力的截斷半徑(r_c)。也就是以中心粒子為球心,r_c為半徑內的所有粒子都是鄰居。skin是比r_c多出來的一部分,所以建立鄰居列表的時候是以r_c+skin為半徑進行尋找的。為什么要這么做呢?這是因為建立鄰居列表是挺費時間的一件事,那么好不容易建立起來那就盡可能多用幾次。而事實上由于分子動力學模擬的時間步長較短,粒子在一個時間步長內走不了多少距離。那么只要鄰居列表中的原子移動距離都小于skin,那么這個列表就可以在下一步中繼續(xù)使用。這就是skin的含義。后面style指的是建立鄰居列表的方法,最常用的就是bin方法。這個方法的具體含義可以查看分子動力學教材。入門只需要知道bin在大多數情況下都是最快的,用它就沒錯了。而r_c在后面的命令中設置。


第五條命令是設置鄰居列表的更新頻率,即neigh_modify命令。lammps提供了很多方式設置鄰居列表的更新頻率。如果鄰居中的粒子已經不在r_c+skin范圍內了,那么鄰居列表就必須要更新了。當需要更新而未更新時就會出現dangerous builds。當算完之后lammps輸出的dangerous builds不為零,那么你就要修改neigh_modify命令了,或者加大skin的值。這條命令使用比較復雜,需要根據模擬經驗和具體的模擬測試確定最優(yōu)值。對于新手,可以直接使用以下設置:

neigh_modify every 1 delay 0 check yes

以上都是計算設置,每個模擬這些設置都大同小異。下面開始就是針對模擬對象的設置。

第六條命令設置原子類型,即atom_style命令。我們模擬的基本構成是原子。不同的對象粒子所具有的屬性是不同的。對于某個對象,粒子不需要的屬性就不需要設置,以節(jié)省內存加速計算。比如研究對象是分子(如水分子),那么構成分子的粒子(style為full)就需要具有以下屬性:粒子編號,粒子所屬分子的編號,粒子類別(atom type),帶電量,x坐標,y坐標,z坐標。如果研究對象是惰性氣體(如氬)或金屬單質(如鐵),那么粒子(style為atomic)的屬性就是:粒子編號,粒子類別,x坐標,y坐標,z坐標。lammps中內置了相當多的原子類型可以查看手冊尋找自己需要的類型。這里要說的是atom style和atom type中文翻譯都一樣,具體是啥意思。簡單來說atom style定義了大類,atom type定義了大類中的小類。比如一個水分子中,所有粒子的atom_style都是full,但是在這個style下,又區(qū)分為氧原子type和氫原子type。一個模擬中可以由多個type,但是只能有一個style。atom_style定義好后就約束了data文件中具體內容。atom_style的語法是

atom_style <style> #舉例:atom_style full

粒子屬性定義完后,就可以定義粒子之間相互作用的方式了。所謂相互作用其實就是粒子之間作用力的計算公式。計算公式是通過勢能來表述的。因為勢能在空間上的梯度就是力,即f=-divu。比如重力勢能對高度求導就是重力。為什么要用勢能表述呢?因為勢能是標量而力是矢量,因此使用標量在構建或者擬合相互作用表達式時就容易的多。而在實際計算中只需要對勢能進行求導就能獲得受力。任何對象的組成粒子之間相互作用都分為非鍵結和鍵結。非鍵結就是不同通過化學鍵的相互作用,鍵結就是通過化學鍵的相互作用。有些對象只有非鍵結相互作用,如液態(tài)氬或單質鐵。有些對象兩種都有,比如水或蛋白質。非鍵結相互作用包括范德華和靜電。在lammps非鍵結相互作用即范德華和靜電通過pair_style和pair_coeff命令設置。pair_style和pair_coeff定義了短程相互作用,即相互作用具有截斷半徑。范德華僅有短程相互作用,而靜電是既有短程又有長程作用,如果體系中有靜電作用,那么定義完pair_style和pair_coeff,還要定義長程靜電求解類型kspace_style。用在分子體系的鍵結相互作用分為鍵(bond),鍵角(angle),二面角(dihedral),離平面(improper)。不同的對象會包含不同的鍵結作用,比如水分子就只包括鍵和鍵角,而分子量稍大的分子會需要全部的鍵結相互作用來描述。這些相互作用構成描述體系的力場。這些力場是從量子力學出發(fā)而得到的經典近似。以下為定義相互作用的命令:

第七條命令設置pair_style。首先定義非鍵結相互作用。不同的力場有不同的pair_style,模擬中要根據對象選擇合適的pair_style。lammps集成了目前幾乎所有的pair_style。模擬時只需要根據對象進行選擇即可。pair_style對于絕大多數模擬都需要定義。不同的pair_style具有不同的語法,具體可查閱手冊。pair_style中一般會定義截斷半徑的大小。以最簡單的包含范德華的pair_style為例,其語法是

pair_style lj/cut <cutoff> #舉例:pair_style lj/cut 9.0

第八條命令設置bond_style。這條命令是可選的,根據不同的對象可定。其設置簡單,沒啥好說的,就是計算兩個粒子通過化學鍵的相互作用,其語法為

bond_style <style> #舉例:bond_style harmonic

第九條命令設置angle_style。這條命令是可選的,根據不同的對象可定。其設置簡單,沒啥好說的,就是計算相鄰三個粒子的相互作用,其語法為

angle_style <style> #舉例:angle_style harmonic

第十條命令設置dihedral_style。這條命令是可選的,根據不同的對象可定。其設置簡單,沒啥好說的,就是計算四個相鄰粒子之間的相互作用,其語法為

dihedral_style <style> #舉例:dihedral_style opls

第十一條命令設置improper_style。這條命令是可選的,根據不同的對象可定。其設置簡單,沒啥好說的,就是計算四個相鄰粒子之間的相互作用,其語法為

improper_style <style> #舉例:improper_style opls

第十二條命令設置special_bonds。這條命令是可選的,當有鍵結作用出現時,設定有鍵結相互作用的粒子之間是否還具有pair相互作用。一般來說,是沒有的,所以這條命令一般寫作:

special_bonds lj/coul 0.0 0.0 0.5

第十三條命令設置kspace_style。這條命令是可選的,當有靜電作用出現時,需要加上。不同的pair_style對應不同的kspace_style,其語法:

kspace_style <style> <value> #舉例:kspace_style pppm 0.0001

以上就是第一部分的內容。對于不同的研究對象力場設置部分會有所不同。如簡單液體只需要pair_style就行了,分子物質一般都會設置,晶體也一般只需要pair_style。如果你的對象有多種物質還可以使用hybrid類的力場命令進行組合設置。

第二部分 定義對象的幾何模型

第十四條命令read_data。這條命令用來讀取粒子的初始位置,鍵結信息等。這些信息寫在一個data文件中。建立data文件稱為lammps的建模過程。lammps的建模是個非常重要的步驟。我們可以在in文件中利用creat_box和create_atoms命令建立。但是這兩條命令只能用來創(chuàng)建晶體結構。對于分子結構需要自己使用第三方工具或者編程建立data文件實現。我習慣于自己寫data文件建立模擬對象的構型。lammps的data文件建立比較麻煩,一般要使用多個工具聯合建模。這里不再贅述,后面會專門講lammps建立data文件的教程。data文件中的信息要根據atom_style和分子的拓撲結構進行建立。以一個由分子構成的模擬對象,也即atom_style是full為例data文件內容為:

lammps description #第一行寫描述信息

116803 atoms #共多少個atoms

70386 bonds #共多少個bonds

41643 angles #共多少個angles

13700 dihedrals #共多少個dihedrals

2550 impropers #共多少個impropers

191 atom types #共多少個atoms types

195 bond types #共多少個bonds types

356 angle types #共多少個angles types

548 dihedral types #共多少個dihedrals types

102 improper types #共多少個angles types

0 97.92 xlo xhi #x方向的上下邊界

0 97.92 ylo yhi #y方向的上下邊界

-15.0 160 zlo zhi #z方向的上下邊界

masses

1 15.9994 # 1號atom type的質量

2 1.008 # 2號atom type的質量

......

atoms

1 1 1 0.0 0.0 0.0 145.0 #粒子編號,粒子所屬分子的編號,粒子類別(atom type),帶電量,x坐標,y坐標,z坐標

2 1 1 0.0 0.0 2.04 147.04

......

bonds

1 2 13825 13826 #bond的編號,bond type,組成bond的第一個,第二個的原子編號

2 3 13826 13827

.......

angles

1 2 13825 13826 13827 #angle的編號,angle type,組成angle的第一個,第二個,第三個的原子編號

2 3 13825 13826 13828

......

dihedrals

1 1 13825 13826 13828 13829 #dihedral的編號,dihedral type,組成dihedral的第一個,第二個,第三個,第四個的原子編號

2 2 13826 13828 13829 13830

......

impropers

1 1 13826 13825 13827 13828 #improper的編號,improper type,組成improper的第一個,第二個,第三個,第四個的原子編號

2 2 13835 13834 13836 13837

......

第三部分 定義對象的物理模型

這部分命令用來指定勢函數中的具體參數是多少。

第十五條命令pair_coeff。根據定義的pair_style和粒子type的個數來決定pair_coeff的內容和條數。具體書寫時可查看手冊中對應的pair_style的例子。如pair_style是lj/cut,粒子type個數為2,則該命令寫為

pair_coeff 1 1 0.1553 3.166 9.0

pair_coeff 1 2 0.1423 3.125 9.0

pair_coeff 2 2 0.2753 3.196 9.0

第十六條命令bond_coeff。根據定義的bond_style和bond type的個數來決定bond_coeff的內容和條數。具體書寫時可查看手冊中對應的bond_style的例子。

第十七條命令angle_coeff。根據定義的angle_style和angle type的個數來決定angle_coeff的內容和條數。具體書寫時可查看手冊中對應的angle_style的例子。

第十八條命令dihedral_coeff。根據定義的dihedral_style和dihedral type的個數來決定dihedral_coeff的內容和條數。具體書寫時可查看手冊中對應的dihedral_style的例子。

第十九條命令improper_coeff。根據定義的improper_style和improper type的個數來決定improper_coeff的內容和條數。具體書寫時可查看手冊中對應的improper_style的例子。

第四部分 定義分組信息

很多情況下我們需要對研究對象的不同部分施加不同的操作。lammps中通過對對象中的原子進行分組實現。

第二十條命令 group。group命令可以采用多種方式進行分組定義

region fix block inf inf inf inf inf 10.0 units box

group fix region fix

group thermos type 1

group a1 id >= 76

group sub id 100:10000:10

group big molecule 123

group polya molecule <> 50 250

group boundary subtract all a2 a3

group boundary union lower upper

group boundary intersect upper flow

group boundary delete

group mine dynamic all region myregion every 100

第五部分 弛豫

在弛豫部分我們首先需要定義基本輸出日志信息,這是為了監(jiān)控我們的模擬是否正確。一般輸出的都是系統量,如溫度、壓強、總動能、總勢能等。所以

第二十一條命令thermo。thermo命令控制每多少步輸出一次。如果不設置lammps默認只在第一步和最后一步輸出日志信息。語法是

thermo <n> #舉例:thermo 1000

第二十二條命令thermo_style。thermo_style控制究竟要輸出哪些全局量。具體能夠輸出哪里量可查閱手冊thermo_style命令描述。如果不設置該條命令,那么lammps就會輸出默認信息。其語法是

thermo_style <style> <args>

第二十三條命令dump。該條命令用來每隔一段時間輸出系統所有原子的信息,這些信息可以是原子的位置,受力,速度,質量,編號等。當你要輸出所有原子的信息時就是用這條命令。將這條命令的輸出文件進行可視化可以初步判斷你的模擬是否正確。以下是最常用的設置

dump 1 all custom 10000 dump.lammpstrj id mol type mass x y z

第二十四條命令minimize。這條命令是可選的。如果你的幾何模型和物理模型都設置正確,那么可以不用這條命令。這條命令執(zhí)行的是能量最小化。為什么要執(zhí)行能量最小化呢?這是因為我們在建模的時候,有時幾何模型和物理模型不是那么完美,經常出現的就是某些原子之間的位置靠得太近,這就會導致系統的能量很大,容易出現運行錯誤。執(zhí)行能量最小化就是將系統的能量調節(jié)至盡可能小。lammps實現這個命令的方式是通過計算原子之間的受力和能量,通過一些優(yōu)化算法,移動原子的位置從而是體系總能量達到盡可能小。所以如果你發(fā)現你的體系中能量超級大,壓強超級大而運行出錯,那么你可以試試能量最小化命令。一般使用下面這些參數設置就能實現你的目的,如果能量最小化后還出錯,那么你就要排查你的幾何模型和物理模型是不是太差了,能量最小化也無能為力了。至于能量最小化采用什么算法,用默認的就行了,也就是說你不用有額外的設置,直接minimize就行了,花里胡哨的沒啥用。

minimize 0.0 1.0e-8 1000 100000

第二十五條命令compute。該命令是lammps中極其重要的一類命令。當你需要計算某個量時你就可以考慮使用compute命令。lammps提供了豐富的compute命令,具體使用時,根據你的需求,查閱手冊看lammps中有沒有實現你需要的功能的compute。compute和dump一樣可以針對某個group的粒子進行計算。一個in文件中可以有多個compute,但每個compute的id不能重復。計算出來的值可以被很多命令使用,如dump,fix等。最常用的就是計算溫度,如

compute myt all temp

第二十六條命令velocity。velocity命令有不少功能,最常用的功能就是根據目標溫度創(chuàng)建粒子的初始速度。最常用的形式為

velocity all create 300.0 4928459 rot yes dist gaussian

第二十七條命令timestep。該命令設定求解牛頓差分方程中時間步的大小。不同體系會要求不同大小的時間步長

。一般來說分子體系會要求時間步長為1fs。簡單液體體系可放寬至5fs。晶體體系可以是2fs。反應力場reaxff體系會要求更小的時間步,0.1-0.5fs。時間步長的選取原子是,要足夠小以準確捕捉體系中最快的動力學過程,又要足夠大以節(jié)省計算時間。通常你的時間步長大小參考文獻設置就行了;虬凑丈厦媾e的例子進行設置,就沒啥問題了。

第二十八條命令fix。fix是lammps中另一類極其重要的命令。其實現功能就是對某個組的粒子施加某種操作。最重要的功能就是用來更新某個組的粒子的速度和位置,因此每個in文件中都必然會有fix命令。我們通常所說的系綜就是通過fix命令實現。如fix nve實現nve系綜,fix nvt實現nvt系綜,fix npt實現npt系綜。關于系綜這里要多說一點。系綜是統計力學的基本概念。所謂系綜,簡單來說就是系統無限個副本的集合。但是分子動力學中的系綜和統計力學上的系綜有聯系又有區(qū)別。比如在lammps中nve系綜就是只直接求解牛頓方程而不對速度和位置做任何干涉的系綜。把fix nve與溫度控制fix聯用就實現nvt系統,把fix nve與溫度控制的fix和壓強控制的fix同時聯用就實現了npt系綜。也可以只使用fix nvt命令直接實現nvt,只使用fix npt命令實現npt系綜。系綜當你的邊界條件是p或f時系統體積是不變的。當你使用s邊界時系統體積是可變得,此時也就沒有什么nvt或nve了。但是需要強調的是npt是通過改變盒子體積來實現控壓的所以在控壓的方向上只能使用p邊界條件。所以基本上你可以忘記系綜的概念。把fix nve或fix nvt或fix npt當做更新原子位置和速度的命令,fix nve直接根據牛二進行更新,fix nvt在更新的時候會對速度進行干涉以達到設定的溫度,fix npt會對速度和位置同時進行干涉以達到設定的溫度和壓強。至于nvt和nve下體積變不變取決于你的邊界條件。同一個系統中可能同時存在fix nve和fix nvt,這是因為有些模擬要求一部分原子的位置和速度直接根據牛二進行更新,而另一部分卻要求在更新原子速度時保持溫度不變。但是不會同時出現fix nve和fix npt,因為fix npt在根據設定的壓強在更新原子位置和速度時,計算的是所有原子對壓強的貢獻。所以fix npt不適合局部控壓。

某個組的粒子只有使用了fix+系綜的命令才會更新位置和速度。因此如果你想保持某組的粒子在模擬過程中固定不動,最簡單的辦法就是不對它施加任何fix命令。通常我們在使用fix時都會使用先來一條fix+系綜,然后調用其他fix實現某些目的,如要讓模擬對象發(fā)生一定變形,就用fix deform,要讓對對象施加一個電場就用fix efield。當你想用實現某種操作時,打開手冊查看fix的一些命令看有沒有實現你功能的fix命令。fix的某些命令也可以用于輸出信息,如fix ave/chunk輸出模擬盒子不同位置處的原子屬性的時間平均值;fix ave/time輸出某些全局量的平均值。

第二十九條命令 run。終于要開始跑了。這個命令最常用的就是你要run多少步。在弛豫階段一般會run 1ns的長度,通常是1000000步。這里可以控制采用什么算法對牛二進行求解。采用默認的就行,也就是說你不用額外設置,直接run就行了。花里胡哨的沒啥用,默認的往往是最好的。其設置是

run 1000000

以上就是弛豫部分。為什么要進行弛豫呢?這是因為系統中原子的初始位置和初始速度都是認為設定的,這是系統的狀態(tài)并不在熱力學平衡狀態(tài)。通常我們希望在正式模擬時系統的起點在熱力學平衡狀態(tài)。只有經過弛豫后系統才會達到熱力學平衡狀態(tài)。在弛豫開始的時候你會看到系統的溫度和壓強會有非常劇烈的變化,這是正常的。經過一段時間后溫度和壓強這些量就會在一個值附近小幅振蕩。那么怎么判斷系統是否達到平衡了呢?這是個挺麻煩的事,不同的體系有不同的標準。通常我們會看系統的溫度,壓強,能量這些值不會在持續(xù)上升或下降就差不多平衡了。但是這個事很難說。所以為了保險起見,通常模擬中會設置冗余的弛豫步數來確保系統達到平衡。這個值一般是1000000步。需要說明的是,系統溫度是根據系統動能算出來的,一般很快就不咋變了。但是壓強是根據受力算出來的振蕩會很大。一般要看溫度或壓強的時間平均值,也就是不同時間步的平均值來判斷系統的溫度或壓強是否達到了設定值。

第三十條命令write_restart。這個命令可以寫出restart文件,以便后續(xù)從這個點接著算。restart文件是一個二進制文件。使用這條命令是個好的習慣,可以讓你省很多事,又不會額外增加計算量。其語法是

write_restart <filename> #舉例write_restart restart.equ

最后一部分 正式模擬

這部分其實才是模擬的核心。不同的研究對象和目的設置和使用的命令都會不同。無法在基本流程中寫,將在后面的專題中介紹。

5. 建模和后處理軟件概述

本節(jié)介紹了包裝 lammps 的商業(yè)和免費軟件,為開發(fā)模型、運行仿真和分析結果提供了一個用戶友好的環(huán)境。

materials design, inc.

materials design, inc. 開發(fā)了medea®,這是一種原子模擬和建模環(huán)境,可提供與lammps一起使用的生產力、模型構建和分析工具。medea®使用流程圖簡化了lammps模擬。生成的lammps模擬很容易在同事之間共享、編輯以供將來重復使用,并且可以由lammps專家進行定制。

使用medea®,可以使用各種晶體系統方法構建原子模型,并使用medea®非晶和熱固性構建器,這是為lammps創(chuàng)建輸入的工具。它們還為預測機械性能提供了經過驗證的方法,包括基于lammps模擬的彈性常數、擴散率、傳輸性能和內聚能。medea®還有助于管理用于模擬有機、無機和金屬系統的力場;材料設計團隊在為有機和金屬系統開發(fā)精確力場方面擁有良好的記錄。

scienomics

scienomics公司開發(fā)了一個lammps接口,作為其材料和工藝模擬 (maps) 平臺的一部分,該平臺允許lammps新手和專家快速有效地創(chuàng)建lammps輸入文件-用于原子模擬和dpd模擬。scienomics為 lammps用戶提供完整的電話和電子郵件支持。maps中的lammps插件還允許用戶為lammps創(chuàng)建輸入文件,對lammps模擬的輸出進行可視化和分析。maps也有許多力場可供選擇(用于原子建模和粗粒度建模)與 lammps。他們還有一個名為amorphous builder和crosslink builder的工具,可用于為lammps創(chuàng)建輸入幾何圖形,以及用于分析reaxff模擬結果的工具和其他強大的功能,用于計算彈性屬性、結構屬性、傳輸屬性等關鍵屬性.

maps是一個綜合性的軟件平臺,結合了構建-模擬-分析的工作流程。這使用戶能夠以高效和無縫的方式對復雜材料系統進行建模,從包括lammps在內的一系列仿真引擎中進行選擇以對過程進行建模,并分析輸出以提取許多屬性,例如機械、熱、傳輸、流變等. 有關所有這些工具和服務的更多詳細信息,請參見scienomics網頁。

xenoview

由sergei shenogin倫斯勒納米技術中心)開發(fā)。該軟件對非商業(yè)用戶免費。 xenoview是基于windows的分子動力學模擬軟件。它的界面為從頭開始構建結構提供了廣泛的工具。此外,您可以從多種格式(例如pdb格式)導入結構。 鍵可以自動定義,力場可以自動分配。xenoview可以導出可在lammps中用于大規(guī)模模擬的數據文件和輸入腳本。有關詳細信息,請參閱xenoview主頁。

atomman

atomman 工具是由nist的lucas hale開發(fā)的開源(免費)軟件。它是一個用于與大規(guī)模原子系統交互的python包。它允許用戶完全從python中準備、運行和分析md模擬。

l 允許對數百萬個原子進行高效和快速的計算,每個原子都有許多自由定義的每個原子屬性。

l 生成包含缺陷的原子系統,例如點缺陷和位錯單極子。

l 用于缺陷分析的內置函數,例如滑移向量和 nye 張量。

l 直接從 python 調用 lammps 并立即檢索結果數據或 lammps 錯誤語句。

l 輕松將系統與 ase 和 pymatgen 的其他 python 原子環(huán)境相互轉換。

l 可以加載基于 cif 晶體結構文件、poscar 和 lammps 原子和轉儲文件的系統。

l 內置單位轉換工具。

l 獨立于平臺:適用于 linux、windows 和 mac。

loos = lightweight object-oriented structure analysis library

loos是一個用于分析分子動力學模擬的軟件包。 它與包無關、開源,可在所有主要的linux發(fā)行版以及osx上運行。它分發(fā)了大約150個預打包的分析工具,范圍從標準任務(軌跡操作、主成分分析等)到新工具(評估模擬收斂、測量膜和膜蛋白的特性)。 此外,它專為快速開發(fā)新的分析工具而設計,尤其是在使用python包裝器時。loos可從https://github.com/grossfieldlab/loos下載

freud = analysis library for post-processing md and mc data

freud是一個python庫,由sharon glotzer的小組(密歇根大學)開發(fā)和支持,它提供了一套簡單、靈活、強大的工具,用于分析從分子動力學或蒙特卡羅模擬中獲得的軌跡。 高性能、并行化 c++ 用于計算標準工具,例如徑向分布函數、相關函數、階參數和聚類,以及包括平均力和扭矩 (pmft)勢和局部環(huán)境匹配在內的原始分析方法。freud庫支持多種輸入格式和輸出numpy數組,支持與許多典型材料科學工作流程的科學python生態(tài)系統集成。

教程:g(r) 分析

文檔:https://freud.readthedocs.io

源代碼:https://github.com/glotzerlab/freud

molecular builders

要模擬分子系統,lammps要求您輸入分子拓撲結構(鍵、角、二面體等列表)以及適合您的模型的力場系數。因此,構建分子系統的任務是一個預處理步驟,并且本身可能是一項復雜的任務。

本節(jié)介紹有助于自動化此過程的工具。

一般來說,封裝可以從鍵拓撲推斷角度、二面角和不正確的相互作用。 他們有生成分子幾何的命令。他們可以從packmol和其他pdb文件生成器生成的文件中讀取坐標。

enhanced monte carlo (emc)

由位于https://verizon.net的't veld (basf)的pieter j.開發(fā)和維護。

增強型monte carlo或簡稱emc為使用compass、charmm、opls、martini、dpd或膠體力場創(chuàng)建和操作粒子模擬的輸入結構提供了環(huán)境。為此,腳本語言管理對其功能的訪問。當前版本通過smiles字符串提供對分子或粗粒度結構的操作,在需要時為選定的力場鍵入這些結構,并應用蒙特卡羅原理構建構象以不重疊原子。emc提供lammps、pdb和xyz格式的輸出端口。可以在https://montecarlo.sourceforge.net上找到適用于linux、macos或windows的編譯版本。

moltemplate

由andrew jewett (ucsb) 開發(fā)和維護。

該工具與 lammps 一起在tools/moltemplate目錄中分發(fā)。有關更多詳細信息,請參見該目錄和moltemplate主頁。

moltemplate專為構建粗粒度生物分子模型而設計。moltemplate可以同時創(chuàng)建:lammps data文件(包含幾何和拓撲)和lammps input腳本(包含力場、修復和組)。與其他轉換工具生成的文件不同,moltemplate允許用戶訪問 lammps中可用的所有力場。用戶可以將分子保存為moltemplate緊湊、可讀的模板文件格式(“.lt”),并與他人共享。分子可以用作更大分子的構建塊。 “罐裝”力場(例如 dreiding、gaff、trappe 和用戶創(chuàng)建)也可以(原則上)以這種格式保存并稍后應用于分子。

分子可以被復制、組合和鏈接在一起以定義新分子。(這些可用于定義更大的分子。支持無限級別的對象組合、嵌套和繼承。)一旦構建,可以自定義單個分子和子單元(原子和鍵,以及子單元可以移動、刪除和替換)。

vmd topotools

由axel kohlmeyer(temple u開發(fā)和維護。

有關詳細信息,請參閱 topotools 主頁。

topotools是一個分子生成器,它利用vmd和tcl的強大功能創(chuàng)建lammps data文件并將它們轉換為其他格式或從其他格式轉換。topotools有兩個組件:一個可以提取和操作拓撲信息的中間件腳本,以及在它之上構建的幾個高級應用程序,例如可以使其能夠讀取/寫入數據文件、復制和合并系統。 與vmd一起,topotools可以從pdb文件、psf文件和原子對距離推斷拓撲結構,使蛋白質溶劑化。

intermol a conversion tool for molecular dynamics simulations

intelmol是用python編寫的,可以在intermol中本地進行 desmond<=>gromacs<=>lammps轉換;amber->x 是通過將amber 轉換為 gromacs,然后使用parmed 轉換為其他程序來執(zhí)行的。amber->charmm 由parmed直接執(zhí)行。

avogadro

作為開源項目開發(fā)和維護。有關詳細信息,請參閱 avogadro 主頁。

avogadro是一款先進的分子編輯器和可視化工具,專為計算化學、分子建模、生物信息學、材料科學和相關領域的跨平臺使用而設計。 它提供靈活的高質量渲染和強大的插件架構。

packmol

作為開源項目開發(fā)和維護。有關詳細信息,請參閱 packmol 主頁。

packmol通過在定義的空間區(qū)域中包裝分子來創(chuàng)建分子動力學模擬的初始點。包裝保證短程排斥相互作用不會破壞模擬。

atomsk

作為開源項目開發(fā)和維護。有關詳細信息,請參閱 atomsk 主頁。

atomsk旨在創(chuàng)建、操作和轉換原子系統。它支持多種文件格式,其中包括 lammps文件格式,還支持vasp、quantum espresso、imd、dl_poly、atomeye cfg格式或xcrysden xsf格式,可以輕松轉換文件以進行從頭計算、經典勢模擬或可視化。 此外,atomsk還可以對原子位置進行一些簡單的變換,如旋轉、變形、插入位錯。

octa 和j-octa

octa是一個開源軟件包,由模擬引擎(分子動力學、流變模擬、自洽場理論、有限元方法等)和用于建模軟物質系統的gui(可視化、簡單分子構建器和分析工具)組成。

octa還為多種模擬器的協作使用提供了環(huán)境,即多物理場和多尺度模擬。

lammps和cognac(md 引擎)文件之間的轉換程序也包括在內。此功能使lammps用戶能夠使用octa gui并結合其他理論(如 scft)運行經典的md模擬。

通過使用商業(yè)版本的j-octa,還可以進行全原子和粗粒度md的復雜分子構建。

molecular simulation design framework (mosdef)

mosdef提供了一組可擴展的python工具,旨在促進使用分子模擬對軟物質系統進行初始化、原子分型和篩選。

automated topology builder (atb) and repository - molecular topology building blocks for organic molecules

atb項目由alan e. mark教授(澳大利亞昆士蘭大學)領導;atb網站上有一份其他貢獻者的列表。

atb以與lammps和其他分子動力學軟件包兼容的格式提供有機分子的拓撲文件。lammps拓撲文件支持使用moltemplate工具(隨lammps分發(fā))構建復雜系統。moltemplate允許結合力場參數和分子模板文件來構建復雜的系統,從而實現類似于在生物分子md模擬包(例如gromos、amber、charmm等)中構建系統拓撲的工作流程。

atb 站點為各種分子(> 20,000并且還在增長)提供拓撲結構。可以使用基于名稱或化學式的內部搜索工具并從結果列表中選擇分子來找到分子。一旦進入感興趣分子的頁面,您可以選擇“molecular dynamics (md) files”并選擇“l(fā)ammps”作為輸出“格式”如果數據庫中不存在分子,可以通過提交新的要處理的分子。

data sites

本節(jié)列出了提供可與lammps一起使用的數據的www站點。例如,原子配置,可用于初始化lammps模擬的模型。或者使用lammps完成的計算存檔,可以瀏覽以獲取物理或模型。

jarvis database for md potential calculations on dft geometries

jarvis for force-fields (jarvis-ff)是一個高通量計算數據庫,用于對具有各種力場/原子間勢的密度泛函理論 (dft)優(yōu)化的幾何結構進行l(wèi)ammps計算。 該項目的目標是通過網絡界面為評估力場提供一個簡單的查找表,并提高數據的可重復性。jarvis-ff是美國國家標準與技術研究院 (nist)材料基因組計劃(mgi)的一部分。

orsi group at queen mary university of london

mario orsi的小組為各種水系統(不同的力場)收集了很好的lammps輸入腳本和數據文件。他們也在為蛋白質和脂質膜系統添加文件。

[ Last edited by 月只藍 on 2022-5-22 at 17:11 ]
回復此樓

» 收錄本帖的淘帖專輯推薦

缺陷計算 科研軟件 分子動力學模擬相關 lammps

» 本帖已獲得的紅花(最新10朵)

» 猜你喜歡

» 本主題相關商家推薦: (我也要在這里推廣)

» 搶金幣啦!回帖就可以得到:

查看全部散金貼

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

luwis

至尊木蟲 (職業(yè)作家)



luyaobao(金幣+5): 謝謝參與
2樓2022-05-22 15:27:52
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

whgregory

木蟲 (正式寫手)



luyaobao(金幣+5): 謝謝參與
求更新,感謝

發(fā)自小木蟲Android客戶端
5樓2022-05-22 17:16:48
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

Lh725

新蟲 (小有名氣)



luyaobao(金幣+5): 謝謝參與
30樓2022-05-30 11:01:03
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

Euphoriagg

新蟲 (小有名氣)



luyaobao(金幣+5): 謝謝參與
31樓2022-06-14 09:53:09
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

haisheng1792

新蟲 (文壇精英)



luyaobao(金幣+5): 謝謝參與
34樓2022-06-24 12:46:05
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

我是Alun

新蟲 (初入文壇)



luyaobao(金幣+5): 謝謝參與
點贊
35樓2022-08-18 10:21:29
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

haibinlucky

木蟲 (小有名氣)



luyaobao(金幣+5): 謝謝參與
厲害!
37樓2022-09-20 11:54:18
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

yptv7361

禁蟲 (著名寫手)

本帖內容被屏蔽

38樓2022-09-23 11:48:54
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

學員0C9LnR

鐵蟲 (小有名氣)



luyaobao(金幣+5): 謝謝參與
點贊!期待專題六!
41樓2022-11-21 15:42:16
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

inaturelz

木蟲 (著名寫手)



luyaobao(金幣+5): 謝謝參與
mark一下,樓主加油

發(fā)自小木蟲Android客戶端
45樓2023-04-20 14:22:41
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖

qiyangtt

鐵蟲 (小有名氣)



luyaobao(金幣+5): 謝謝參與
樓主有興趣做代算嗎,+微tjs7188私聊
49樓2023-11-09 14:03:20
已閱   回復此樓   關注TA 給TA發(fā)消息 送TA紅花 TA的回帖
簡單回復
2022-05-22 15:50   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-22 15:58   回復  
luyaobao(金幣+5): 謝謝參與
tzynew6樓
2022-05-22 18:37   回復  
luyaobao(金幣+5): 謝謝參與
i 發(fā)自小木蟲Android客戶端
gter_wang7樓
2022-05-22 19:13   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-22 19:31   回復  
luyaobao(金幣+5): 謝謝參與
yes 發(fā)自小木蟲Android客戶端
nono20099樓
2022-05-22 21:39   回復  
luyaobao(金幣+5): 謝謝參與
?
gter_wang10樓
2022-05-22 22:01   回復  
tigerwood211樓
2022-05-22 23:12   回復  
luyaobao(金幣+5): 謝謝參與
ok
jiaoxg12樓
2022-05-23 09:14   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-23 16:45   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-23 17:13   回復  
luyaobao(金幣+5): 謝謝參與
5 發(fā)自小木蟲IOS客戶端
kanghuijun15樓
2022-05-23 18:41   回復  
luyaobao(金幣+5): 謝謝參與
dsctg16樓
2022-05-24 08:26   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-24 09:10   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-24 16:43   回復  
luyaobao(金幣+5): 謝謝參與
假大空19樓
2022-05-24 19:23   回復  
luyaobao(金幣+5): 謝謝參與
tigerwood220樓
2022-05-25 06:36   回復  
ok
2022-05-25 08:52   回復  
luyaobao(金幣+5): 謝謝參與
1
gdflly0822樓
2022-05-25 21:48   回復  
luyaobao(金幣+5): 謝謝參與
2022-05-26 09:17   回復  
luyaobao(金幣+5): 謝謝參與
kanghuijun24樓
2022-05-26 13:57   回復  
2022-05-26 19:48   回復  
luyaobao(金幣+5): 謝謝參與
發(fā)自小木蟲Android客戶端
2022-05-27 09:03   回復  
luyaobao(金幣+5): 謝謝參與
liang888727樓
2022-05-27 14:41   回復  
luyaobao(金幣+5): 謝謝參與
祝福。
myttiger28樓
2022-05-28 10:03   回復  
luyaobao(金幣+5): 謝謝參與
twenty4929樓
2022-05-28 19:56   回復  
luyaobao(金幣+5): 謝謝參與
110001180532樓
2022-06-24 00:02   回復  
luyaobao(金幣+5): 謝謝參與
! 發(fā)自小木蟲Android客戶端
2022-06-24 08:07   回復  
luyaobao(金幣+5): 謝謝參與
發(fā)自小木蟲Android客戶端
XG-WUST36樓
2022-08-18 22:11   回復  
2022-10-17 17:04   回復  
luyaobao(金幣+5): 謝謝參與
XG-WUST40樓
2022-10-17 23:08   回復  
XG-WUST42樓
2022-11-21 22:22   回復  
2023-02-10 16:17   回復  
XG-WUST44樓
2023-02-10 21:54   回復  
luyaobao(金幣+5): 謝謝參與
q 發(fā)自小木蟲IOS客戶端
一剪梅946樓
2023-09-07 14:49   回復  
luyaobao(金幣+5): 謝謝參與
送紅花一朵
f
冰墩墩+48樓
2023-09-09 17:54   回復  
luyaobao(金幣+5): 謝謝參與
2024-02-18 11:53   回復  
luyaobao(金幣+5): 謝謝參與
發(fā)自小木蟲Android客戶端
相關版塊跳轉 我要訂閱樓主 luyaobao 的主題更新
提示: 如果您在30分鐘內回復過其他散金貼,則可能無法領取此貼金幣
普通表情 高級回復 (可上傳附件)
最具人氣熱帖推薦 [查看全部] 作者 回/看 最后發(fā)表
[考研] 學碩材料275調劑 +6 路三三 2026-03-03 6/300 2026-03-04 09:22 by zhyzzh
[考研] 304分材料專碩求調劑 +5 qiuzhigril 2026-03-03 7/350 2026-03-04 08:49 by EBSD
[考研] 0703化學調劑 +3 G212 2026-03-03 4/200 2026-03-04 08:41 by houyaoxu
[考研] 歡迎采礦、地質、巖土、計算機、人工智能等專業(yè)的同學報考 +8 pin8023 2026-02-28 11/550 2026-03-04 00:12 by 不吃魚肉
[考研] 接收調劑 +15 津萌津萌 2026-03-02 23/1150 2026-03-03 23:39 by 夢—-
[考研] 0703化學 學碩 理工科均可 不區(qū)分研究方向 總分279求調劑 +5 1一11 2026-03-03 5/250 2026-03-03 23:14 by zhukairuo
[考研] 292求調劑 +3 sgbl 2026-03-03 3/150 2026-03-03 17:30 by barlinike
[考研] 278求調劑 +3 滿天星11_22 2026-03-02 3/150 2026-03-03 13:51 by Iveryant
[考研] 291求調劑 +3 MuoLuo1312 2026-03-02 6/300 2026-03-03 12:13 by 熱情沙漠
[考博] 26申博 +4 north, 2026-02-28 4/200 2026-03-03 10:12 by LOYUBUAA
[考研] 清華大學 材料與化工 353分求調劑 +5 awaystay 2026-03-02 6/300 2026-03-03 09:03 by houyaoxu
[考研] 11408,學碩276求調劑 +3 崔wj 2026-03-02 5/250 2026-03-02 22:14 by 255671
[考研] 285求調劑 +9 滿頭大汗的學生 2026-02-28 9/450 2026-03-02 20:29 by hypershenger
[考研] 一志愿華南理工大學材料與化工326分,求調劑 +3 wujinrui1 2026-02-28 3/150 2026-03-02 16:36 by chuocheng
[考研] 0856化工專碩求調劑 +15 董boxing 2026-03-01 15/750 2026-03-02 15:06 by 晃晃不許晃
[考研] 化工專碩342,一志愿大連理工大學,求調劑 +6 kyf化工 2026-02-28 7/350 2026-03-02 10:56 by 無際的草原
[考研] 材料類求調劑 +11 wana_kiko 2026-02-28 14/700 2026-03-02 08:46 by 聰明的大松鼠
[考研] 311求調劑 +6 亭亭亭01 2026-03-01 6/300 2026-03-01 15:41 by 324616
[論文投稿] Optics letters投稿被拒求助 30+3 luckyry 2026-02-26 4/200 2026-03-01 09:06 by babero
[考研] 304求調劑 +3 52hz~~ 2026-02-28 5/250 2026-03-01 00:00 by 52hz~~
信息提示
請?zhí)钐幚硪庖?/div>
在线 制服 中文字幕 日韩| 欧美日韩亚洲tv不卡久久| 2026天天操天天干| 国产激情在线观看一区二区三区| 神马不卡视频在线视频| 天天操天天日天天插天天舔| 亚洲永远av在线播放| 亚洲综合在线视频在线播放| 神马午夜久久电影网| 亚洲高清一区二区三区久久| 2018中文字字幕人妻| 亚洲欧美激情国产综合久久久| 成人做爰av在线观看网站| 99免费观看在线视频| 抽插小穴啊啊啊视频| 黄色片免费网站在线| 91超精品碰国产在线观看| 国产精品成人免费电影| 亚洲成人,国产精品| 久久久久久免费观看av| 中文字幕在线字幕乱码怎么设置 | 日本黄色一级电影网址| 抽插小穴啊啊啊视频| 妈妈的朋友中字在线免费观看| 大陆中文字幕视频在线| 久久久久久高清一区| 手机视频在线观看一区| 国产高清自拍偷拍在线| 国产成人av在线你懂得| 麻豆午夜激情在线观看| 亚洲天堂色综合久久| 2020年亚洲男人天堂网| 在线能看视频你懂的| a级黄片免费观看| 2020国产激情视频在线观看| 色网站在线观看免费| 久久久国产精品免费视频网| 中文字幕 首页 人妻| 91精品国产91久久久久久密臀| 中文字幕精品人妻久久久久 | 日本亚洲午夜福利一区二区三区| 国产农村乱子伦精精品视频| 亚成区一区二区人妻熟女| av在线播放观看h| 自拍偷自拍亚洲精品10p| 青青青青午夜手机国产视频| 91大神福利视频网| 欧美日韩综合精品无人区| 污视频在线观看地址| 亚洲人精品午夜射精日韩| 亚洲成人激情在线综合| 久久99嫩草99久久精品| 男人的天堂aⅴ在线| 久久精品国产亚洲av热软件| 日韩最近中文在线观看| 免费在线小视频你懂的| 日本少妇人妻凌辱在线| 亚洲一区二区三区四区入口| 亚洲国产精品自产拍在线观看| 国产三级自拍视频在线观看网站| 不卡一区二区视频在线| 亚洲第一区av中文字幕| 欧美日韩不卡视频合集| 一区二区三区四区久久久久韩日 | 婷婷色九月综合激情丁香| 熟女人妻少妇一区二区| 亚洲春色av中文字幕| 精品国产av虐杀两警花| 神马午夜久久电影网| 亚洲国产综合久久精品| 天天操天天搞天天操| 少妇熟女天堂网av| 亚洲欧美精品日韩偷拍| 欧美日本国产一区二区| 加勒比不卡在线视频| 亚洲精品久久久人妻| 麻豆国产精品777777在| av资源中文字幕在线观看| 天天早上头和脸出汗是怎么办| 青青操久久综合激情| 亚洲欧美精品日韩偷拍| 中文字幕 一区二区在线观看| 亚洲一区视频中文字幕在线播放| 黑人3p日本女优中出| 亚洲色图日韩在线视频观看| 天海翼亚洲一区在线观看| 亚洲综合一区二区三区四区| 天天色 天天操 天天好逼| 青青操久久综合激情| 色丁香久久激情综合网| 亚洲成人偷拍自拍在线| 婷婷色九月综合激情丁香| 中文字幕 中文字幕 亚洲| 无人区一码二码三码区别在哪 | 亚洲男人的天堂最新网址| 放荡人妻极品少妇全集| 交换的一天中文字幕在线视频| 久久99精品久久久久久三级| 久久视频 在线播放| 国产福利一区二区三区在线观看| 美女网站视频久久精品| 新香蕉视频香蕉视频2| 好看的日本中文字幕在线观看二区| 公侵犯人妻中文字幕巨| 高潮喷水在线视频观看| 青娱乐免费最新视频| 亚洲一级熟妇丰满的女人| 91超精品碰国产在线观看| 91人妻人人做人人爽高清| 手机看片福利一区二区三区四区| 夜夜操夜夜爱夜夜摸| 顶级欧美色妇xxxx| 人妻激情偷乱一区二区三区av| 美女福利视频一区二区三区四区| 亚洲国产精品久久久久久无码| 亚洲av激情综合网| 青青在线视频看看| 999久久久人妻精品一区 | 丰满人妻被猛烈进入中文字幕| 大香蕉尹人在线最新| 桃色成人开心激情网| 国产精品蝌蚪自拍视频| 黑人爆操女人免费视频| 东京热男人的天堂视频| tushy一区二区三区视频| 亚洲成人三级黄色片| 精品国产无乱码一区二区三区| 亚洲人妻系列在线视频| 日韩欧美中文字幕老司机三分钟| 天天干天天色综合久久| 色视频免费观看网址| 国产精品剧情av在线播放| 日本清纯中文字幕版| 国内精品一区二区2021在线| 无码人妻丰满熟妇区五路| 东京热日韩av影片| 9久re热视频在线精品| 啊不行啊操逼好爽大鸡吧视频| 2019年中文字幕在线播放视频| 美女网站视频久久精品| 超碰在线观看97资源| 美女福利视频一区二区三区四区| 韩国一级片最火爆中文字幕| 二十四小时日本高清在线观看 | 美女把逼扒开让男人桶| 色欲AV蜜桃一区二区三| 国产做A爱免费视频在线观看| 天堂网免费在线电影| 欧美成人性生活视频播放| 中文字幕欧美人妻在线.| 青青操91美女国产| 99热99这里免费的精品| 欧美成人久久久桃色aa| 夜色17s精品人妻熟女av| 国内精品一区二区2021在线 | 最新激情中文字幕视频| 午夜情色一区二区三区| 38av一区二区三区| 91亚洲国产成人久久精品| 日韩一级视频一区二区三区| 青青青在线观看国产| 92午夜免费福利视频www| 亚洲黄色免费在线观看网站| 国产精品久久久久久成人久| 50熟妇一区二区三区| 黑人巨大精品一区二区在线 | av里面的动作是真进去吗| 中文字幕观看中文字幕免费| 色狠狠色综合久久久绯色| 久草视频在线看免费| 69国产精品成人aaaaa片| 国产三级自拍视频在线观看网站| 国产在线观看一区二区三区四区| 亚洲国产精品青青草| 河北全程露脸对白自拍| 青青青在线视频免费播放| 夏目彩春av在线看| jiee日本美女视频网站| 国产精品黄色片大全| 亚洲综合熟女乱中文| 超peng视频在线免费播放97| 波多野结衣在线一区别| 日韩欧美一区二区三区免费看| 中文字幕熟女人妻丝袜丝在线| 久久99精品热在线观看| 麻豆白洁少妇在线播放| 猫咪亚洲中文在线中文字幕| 黄色大片一级老太太操逼| 不卡一二三区别视频| 黑人侵犯人妻森泽佳奈| 久久免费视频ww一区| 久久久久久a女人处女| 91九色国产在线视频| 色欲AV亚洲AV无码精品| 日本在线免费观看国产精品| 日本一道中文字幕99| 呻吟求饶的人妻中文字幕| 夜色17s精品人妻熟女av| 污网址在线观看视频| 亚洲av手机免费在线| 福利在线国产小视频| 黑人3p日本女优中出| 18禁网站在线点击观看| 国产伦理二区三区在干嘛呢| av福利免费体验观看| 国模伊人久久精品一区二区三区| 久久久久久久精品乱码| av网页免费在线观看| 一区二区三区四区久久久久韩日| 精品免费一区二区三区四区视频| 操死你美女在线视频| 欧美强奸视频在线观看| 日本一区二区三区的资源| 交换的一天中文字幕在线视频 | 久久精品久久久久观看99水蜜桃| 制服丝袜 中文字幕 日韩| 日韩一级欧美一级片| 日韩人妻一区二区三区在线观看| 欧美强奸视频在线观看| 人妻少妇精品二三区| 狠狠操av一区二区三区| 伊人精品成人综合网| 一区二区在线观看视频观看| 国产精美视频精品视频精品| 制服丝袜中文字幕熟女人妻| 国产igao激情在线视频入口| 国产高清自拍偷拍在线| 久久人妻人人草人人爽| 绿巨人浩克在线视频观看 | 中文字幕熟女人妻丝袜丝在线| 日韩人妻中文字幕区| 92午夜免费福利视频www| 大尺度av毛片在线网址| 日韩久久九九精品视频| 91精品久久久久久久99蜜月| 国产成人在线观看视频播放| 伦理在线观看未删减中文字幕 | 午夜精品视频免费观看| 欧美日本在线免费视频| 少妇熟女天堂网av| 公侵犯人妻中文字幕巨| 天天曰天天摸天天爽| 国产激情免费在线视频| 亚洲美女午夜激情视频在线观看| 福利一二三在线视频观看| 国内自拍第一区二区三区| 天天爽天天操天天插| 老牛影视在线一区二区三区| 上床啪啪啪免费视频| 国产伦理二区三区在干嘛呢| 91香蕉国产亚洲一二三区| 国长拍拍视频免费孕妇| 亚洲码av一区二区三区| 凹凸视频一区二区在线观看| 亚洲精品激情视频在线观看| 360偷拍蜜桃臀69式| 日本少妇丰满大bbb的小乳沟| 少妇被中出一区二区| 99国产精品久久99久久久| 青娱乐免费最新视频| 每日更新日韩欧美在线| 麻豆出品视频在线观看| 欧美在线观看一区二区不卡| 2019年中文字幕在线播放视频| 男女插鸡巴视频软件| 夜夜人人干人人爱人人操| 青青青在线观看国产| 人妻激情综合久久久久蜜桃 | 夜夜操天天干夜夜操| 69视频在线精品国自产拍| 最新国产精品拍在线观看| 五月的婷婷综合视频| 欧美一区二区三区视频看| 日本熟妇乱妇熟色视频| 得得爱在线视频观看| 精久久久久久久久久久久| 亚洲午夜熟女在线观看| 免费看日韩黄视频在线观看| 色视频免费观看网址| 神马不卡视频在线视频| 最新福利二区三区视频| 蜜桃tv一区二区三区| 久久久久九九九九九12| 九九热在线精品播放| 国产,亚洲,欧美综合| 东京热日韩av在线| 婷婷六月天在线视频| 男人和女人的逼视频| 中日韩又粗又硬又大精品| 天天日夜夜操人人爽| 日韩男女视频网站在线观看| 人妻激情综合久久久久蜜桃 | 女人扒开逼让男人操| 38av一区二区三区| 天天天天天天天天干夜夜| 天天干天天日天天弄| 天堂av国产av伦理av| 人妻人妻在线视频网站| 天天天天天天天天日日日| 亚洲国产美女主播在线观看| 一区二区三区av免费天天看| 成人精品影视一区二区| 国产精品剧情av在线播放| 91精品综合久久久久久五月天| 91偷拍被偷拍在线播放| 岛国av成人午夜高清| 91色哟哟视频在线观看| 国产精美视频精品视频精品 | 日本东京热视频欧美视频| 亚洲成人 国产精品| 99女福利女女视频在线播放| 羞羞漫画无限免费观看秋蝉| 精产国品一二三产品区别97| 大鸡扒操大逼大片免费关看| 亚洲成人自拍图片网站 | 婷婷色九月综合激情丁香| 亚洲熟妇在线视频观看| 天天操天天干加勒比久久| 国产欧美福利在线观看| 亚洲自拍偷拍一区二区中文字幕 | 亚洲一区二区在线激情| 国产精品igao为爱寻找激情| 亚洲欧美一级特黄大片| 天堂av国产av伦理av| 2020国产成人精品视频| 中文字幕av特黄毛片| 欧美一区日韩二区三区四区| 日韩男女视频网站在线观看| 欧美视频免费观看777| 偷拍欧美日韩另类图片| 亚洲韩精品一区二区三区| 午夜夫妻性生活视频| 99久久久久久久久久久久久| 日韩一区二区在线播放观看| 综合激情网,激情五月| 成年人免费黄色av| 精产国品一二三产品区别91| 伊人网在线观看 视频一区| 69精品人妻久久久久久久久久久| 欧美最新一区二区三区| 亚洲欧美国产一本综合首页| 日本老熟老熟妇七十路| 亚洲永远av在线播放| 男人av一区二区三区| 日本老熟妇av老熟妇| 老司机免费视频福利0| 久久99国产中文丝袜| 大片a免费观看在线视频观看| 男生用大肌巴操美女骚穴| 亚洲一区二区偷拍女厕所| av在线免费在线观看| 九九热在线精品播放| av福利免费体验观看| 在线免费观看欧美小视频| 大尺度av毛片在线网址| 黄色av网址在线播放| 亚洲制服丝袜在线看| 日韩最近中文在线观看| 亚洲午夜高清在线观看| 国产成人情侣av在线| 河北全程露脸对白自拍| 熟女人妻aⅴ一区二区三| 女人的天堂av在线网| 92麻豆一区二区三区| 亚洲欧美小说中文字幕| 老司机免费视频福利0| 黄色av日韩在线观看| 天天插天天干天天狠| 99久久国产精品免费热| 91污污在线观看视频| 久久久精品人妻无码专区不卡 | 国际精品熟女一区二区| 开心激情五月天作爱片| 杜达雄啪啪毛片视频| 91精品麻豆91夜夜骚| 五月激情婷婷四射基地| 久久精品国产亚洲av清纯| 美利坚合众国av天堂| 松本菜奈实最新av在线| 亚洲中文字幕最新地址| 美女精品久久久久久久久| 91大神在线免费观看视频| 亚洲 综合 欧美 一区| 中文字幕欧美一区二区视频| 欧美亚洲另类精品第一页| 最近在线中文字幕免费| 国产美女视频带a∨黄色片| 天天躁狠狠躁狠狠躁性色| 日本黄页在线观看视频| 亚洲综合熟女乱中文| 国产精品视频网站污污污 | 亚州av嫩草av极品在线观看| 人妻视频网站快射视频网站| 91大神福利视频网| 91精品综合久久久久久五月天| 九九热视频1这里只有精品| 亚洲另类欧美综合久久| 亚洲AV无码一二三四区在线播放| 天天日夜夜操人人爽| 日韩激情文学在线视频| 97视频538在线观看| 国产视频1区2区3区| 青娱乐免费视频一二三| 天天操天天干天天谢| 久久sm人妻中出精品一区二区| 久久久久久久精品乱码| 精品av天堂毛片久久久| 日本五六十路熟女视频| 另类欧美激情校园春色| av网页免费在线观看| 性感美女人妻久久久| 亚洲精品1卡2卡3卡| 68福利精品在线视频| 一区二区三区午夜福利在线| 九一精品人妻一区二区三区| 自拍偷自拍亚洲精品10p| 西野翔人妻中文字幕中字在| 久久精品国产亚洲av清纯| 天天日夜夜操人人爽| 高潮喷水在线视频观看| 精品美女洗澡一区二区| 亚洲一区二区中文字幕久久 | 福利美女视频在线观看| 中文字幕熟女人妻一区| 五月激情婷婷四射基地| 91福利高清在线播放| 亚洲欧美激情国产综合久久久| 人妻在线中文视频视频| 亚洲熟女一区二区六区| 9420高清视频在线观看国语版| 大乳丰满人妻中文字幕韩国hd| 亚洲一区二区精品在线播放| 天天插天天干天天狠| 人人妻人人狠人人爽| 国产毛片特级Av片| 日本福利视频网站导航| 亚洲精品激情视频在线观看| 免费的啪啪视频软件| 欧洲成熟女人色惰片| 欧美视频亚洲视频在线| 久久久久久高清一区| 日韩男女视频网站在线观看| 熟女俱乐部jukujoclub| 天天弄天天草天天日天天| 女生裸体视频免费网站| 日本男女免费福利视频| 大奶熟妇激情操逼逼| 91精品夜夜夜一区二区| 99久久99九九九99九| 欧美在线观看一区二区不卡| 国产在线观看av一区| 欧美大胆a级视频秒播| 韩国一级片最火爆中文字幕| 最近最新欧美日韩精品| 人妻激情综合久久久久蜜桃 | 黑人3p日本女优中出| 国产毛片特级Av片| 亚洲美女黄色福利视频网站大全| 三级欧美日韩一区二区三区| 92午夜免费福利视频www| 中文字幕综合网91| 免费在线观看亚洲福利| 美女精品久久久久久久久| 在线播放 日韩 av| 大香蕉在线欧美在线视频| 日本美女爱爱视频网站| 亚洲人成小说网站色| 国产免费久久精品99re丫丫| 青青青在线视频免费播放| 欧美一级aaaaaaa片| 91色哟哟视频在线观看| 成人大片男人的天堂| av激情四射五月婷婷| 九九热精品视频在线播放| 99精品久久一区二区| 97精品久久久久久无码人妻| 夜夜躁av麻豆男| 18岁禁一二三区免费体验| av资源中文字幕在线观看| 亚洲美女a级黄色在线播放| 红桃视频国产av在线| 一区二区三区高清视频3| 蜜桃臀av在线一区二区| av 一区二区三区 熟女| 大成色亚洲一二三区| 日韩av熟妇在线观看| 久久sm人妻中出精品一区二区| 美女张开腿给男人桶爽的软件| 国产福利一区二区三区在线观看 | 国产精品久久久久精品三级18| 亚洲免费午夜污福利| 韩日一级人添人人澡人人妻精品| 亚洲自拍偷拍av在线| 乌克兰美女操逼高清内射视频| 亚洲熟女一区二区六区| 伊人情人成综合视频| 日本少妇丰满大bbb的小乳沟| 91大神福利视频网| 38av一区二区三区| 中文字幕免费啪啪啪| 免费看超污视频在线观看| 亚洲一区二区三区无码在线| 色欲AV亚洲AV无码精品| 亚洲a级视频在线播放| 天天在线播放日韩av| 五十岁熟妇高潮喷水| 欧美亚洲愉拍一区二区三区| 久久国产精品久精国产爱 | 日韩黄色在线观看网站上 | 日本亚洲精品视频在线观看| 国产av嗯嗯啊啊av| 在线观看免费啪啪啪| 得得爱在线视频观看| 久久久久国产精品二区| 亚洲男人天堂最新网址大全| 中文字幕人妻精品精品| 亚洲三级综合在线观看| 熟女人妻aⅴ一区二区三| 色视频在线播放免费观看| 亚洲综合第一区二区| 午夜精品久久秘?18免费观看| 老熟妇一区二区三区v∧88| av里面的动作是真进去吗| 熟女人妻精品视频一区| 国产精品久久人人添| tushy一区二区三区视频| 日本有码精品一区二区三区| 天天躁狠狠躁狠狠躁性色| 亚洲va999天堂va| 黄色片免费网站在线| 99久久国产精品免费消防器材| 亚洲图片另类综合小说| 亚洲成人激情在线综合| 日韩av熟妇在线观看| 亚洲欧美小说中文字幕| 97人妻av人人澡人人爽| 大成色亚洲一二三区| 蜜桃臀av在线一区二区| 插鸡视频免费网站在线播放| 天天摸天天干夜夜操| 在线中文字幕人妻av| 日韩人妻精品久久久久| 网友自拍第一页99热| 5566熟女人妻人妻| 5d蜜桃臀女无痕裸感| 国内精品一区二区2021在线| 国产白丝一区二区三区av| 欧美成人红桃视频在线观看| v天堂国产精品久久| 亚洲a级视频在线播放| 夜夜人人干人人爱人人操| 国产精品免费看一区二区三区| 真人一进一出抽搐大尺度视频| 日韩成人免费观看电影| 国产最新av在线免费观看| av激情四射五月婷婷| 免费中文三级在线观看| av一区二区三区四区五区在线| 国产免费久久精品99re丫丫| 国产三级自拍视频在线观看网站| 欧美精品一区二区三区观看| 最新日韩中文字幕免费在线观看| 日日躁夜夜躁狠狠操| 亚洲午夜精品一级毛片app| 18禁男女啪啪啪无遮挡| 美女一区二区四区六区八区| 55夜色66夜色亚洲精品| 中文字幕久久久国产| 中文字幕一区二区三区久久久| 成人午夜麻豆大胆视频| 不卡一区二区视频在线| 国产精品成人免费电影| 日本不卡 中文字幕| 新香蕉视频香蕉视频2| 18禁男女啪啪啪无遮挡| 久久久久久久精品乱码| 乱子伦国产一区二区三区 | 神马不卡视频在线视频| 神马不卡视频在线视频| 91麻豆精品国产在线| 亚洲成人偷拍自拍在线| 亚洲av综合av一去二区三区| 中文字幕欧美一区二区视频| 亚洲一区二区三区无码在线| 欧美大胆a级视频秒播| 91超碰九色porny| 久久久精品人妻无码专区不卡 | 51vv精品视频在线观看| 久久内射天天玩天天懂色| 182tv精品免费在线观看| 男人的天堂在线2025| 久久久久九九九九九12| 国产成人91色精品免费看片| 国产高清视频www夜色资源| 欧美区日本区国产区| 超碰在线pro中文字幕| 人妻少妇视频系列视频在线| 一区二区在线观看视频观看| 国内销魂老女人老泬| 国产 亚洲 欧美 自拍| 久草视频在线看免费| 亚洲一区二区三区国产精品电影 | 五月天天堂视频在线| 亭亭五月天在线观看| 欧美亚洲国产一区二区| 在宿舍强奷两个清纯校花| 极品风骚人妻3p视频| 一区二区三区国产在线成人av| 成人午夜麻豆大胆视频| 日本不卡视频一二三区| 美女把逼扒开让男人桶| 亚洲午夜国产末满十八岁勿进网站| 91超碰九色porny| 福利视频免费在线播放| 青青青免费手机视频在线观看| 国产女人18毛片水真多精选| 日本老女人日比视频| 啪啪啪网站免费在线看| 国产资源网站在线播放| 国产又粗又长又大视频| 干逼又爽又黄又免费的视频| 自拍偷自拍亚洲精品10p| 日本a级2020在线观看| 河北全程露脸对白自拍| 日本亚洲精品视频在线观看| 中文字幕国产一区在线视频| 大尺度久久久久久久| 亚洲成年人精品国产| 美女张开腿给男人桶爽的软件| 日韩成人免费观看电影| 视频在线+欧美十亚洲曰本| 在线视频国产精品欧美| 小妹妹爱大棒棒免费观看视频| 日韩久久不卡免费视频| 欧美亚洲国产一区二区| 天天操天天舔天天射天天日天天干| 自拍丝袜国产欧美日韩| 亚洲制服丝袜在线看| 国产不卡免费在线观看| 中文字幕在线字幕乱码怎么设置 | 欧美色视频网址大全| 亚洲熟女人妻自拍在线视频| 99热99这里免费的精品| 人人妻人人爽人人摸| 一区二区在线观看视频网站| 天天日夜夜操人人爽| 日本欧美高清在线观看视频| 久久sm人妻中出精品一区二区| 亚洲精品色图1234| 亚洲欧美激情久久久| 日本一区二区三区调教性奴视频| 青青青在线视频免费播放| 在线观看中文字幕视频成人| 18禁网站在线点击观看| 91精品麻豆91夜夜骚| 欧美一区二区三区视频看| 68福利精品在线视频| 欧美久久一区二区伊人| 九九九九九久久久国产| 不卡一二三区别视频| 神马不卡视频在线视频| 国际日韩日韩日韩日韩日韩| 操人妻人妻天天爽天天偷| 内地精品毛片在线观看| 天天色 天天操 天天好逼| 顶级欧美色妇xxxx| 欧美日韩精品aaa| 99精品视频在线在线观看| 亚洲av三级电影在线观看| 一区二区三区四区久久久久韩日| 亚洲蜜桃久久久久久| 亚洲av综合av一去二区三区| av大尺度一区二区三区| 精品久久久久久久久久久久久| 国产,亚洲,欧美综合| alisontyler和黑人| 熟女一区二区三区综合| 亚洲国内精品久久久久久久 | 懂色av之国产精品| 美国男的操女孩的小嫩逼| 亚洲第一成年偷拍视频| 强乱人妻中文字幕日本| 天天爱天天日天天爽| 国产肥胖熟女又色又爽免费视频| 另类欧美激情校园春色| 最近最新欧美日韩精品 | 日韩激情亚洲国产欧美另类激情 | 五月激情婷婷四射基地| 国产人妻777人伦精品hd超碰| 国产成人在线观看hd| 最新激情中文字幕视频| 欧美成人屋影院在线视频观看| 精品人妻人人做人人爽| 亚洲天堂色综合久久| avgo成人短视频| 美女扒开逼逼给你看| 亚洲精品久久久人妻| 精品欧美乱码久久久| 汤姆提醒30秒中转进站口| 三级欧美日韩一区二区三区| 天天操天天干天天舔天天| 国产一区两区三区福利小视频| 七色福利视频在线观看| 女人高潮潮呻吟喷水网站| 亚洲蜜桃久久久久久| 99热99这里免费的精品| 美女把腿张开给男的捅| 99久久国产精品免费热| 五月婷婷激情视频网| 夜色福利视频免费观看| 杜达雄啪啪毛片视频| 午夜国产成人精品视频观看| 偷拍熟女大胆免费视频| 成人超碰一区二区三区| 日本欧美高清在线观看视频| 中文字幕 一区二区在线观看| 亚洲国产精品自产拍在线观看| 男人用大鸡巴狂操女人肉穴| 91九色pony蝌蚪| 免费啪啪啪网站在线观看| 亚洲gay视频在线观看| 亚洲一区二区三区四区入口| 国产激情一区二区视频| 顶级欧美色妇xxxx| ysl蜜桃色7425| 日本人妻熟妇丰满成熟HD系列| 三级欧美日韩一区二区三区| 午夜美女福利视频在线| 亚洲字幕一区二区夜色av| 91精品一区一区三区| 99精品视频在线在线观看| 亚洲永远av在线播放| 97视频538在线观看| 国产大桥未久一区二区| 上床啪啪啪免费视频| 插鸡视频免费网站在线播放| 亚洲自拍偷拍av在线| 国产在线小视频一区二区| 中文字幕人妻一区二区视频系列| 一区二区三区高清视频3| 日本免费人爱做视频在线观看不卡 | 中文字幕av特黄毛片| 91亚洲最新蜜桃在线| 欧美极品少妇高潮喷水| 青青草一个释放的网站| 国产精品成人免费电影| av日韩视频在线观看| 99热在线只有的精品| 亚洲成年人精品国产| 91青青青国产免费高清| 91九色pony蝌蚪| 亚洲一区二区三区国产精品电影| 国产igao激情在线视频入口| 91人妻人人做人人爽高清| 国际精品熟女一区二区| 91香蕉国产亚洲一二三区| 中文字幕人妻一区二区视频系列 | 91精产国品一二三产区区别网站| 91久久久久久最新网站| 99国产精品久久99久久久| 亚洲永远av在线播放| 欧美性受黑人猛交裸体视频 | 精久久久久久久久久久久 | 50熟妇一区二区三区| 伊人精品成人综合网| 美国男的操女孩的小嫩逼| 伊人精品成人综合网| 中文字幕精品人妻久久久久| 美女福利视频一区二区三区四区| 在线免费观看视频18| 欧美性感美女热舞视频| av 资源在线播放| 在线免费观看a视频免费| 二十四小时日本高清在线观看| 日本五六十路熟女视频| 亚洲中文字幕最新地址| av中文字幕国产精品| 91青青青国产免费高清| 91日本精产品一区二区三区| 人妻被强av系列一区二区| 中文字幕熟女人妻一区| 欧美精品熟妇免费在线| 四季av人妻一区二区三区| 免费中文三级在线观看| 欧美最新一区二区三区| 99 re国产精品| 国产自拍偷拍在线精品| 4438x亚洲最大的成人| 99久久人人爽亚洲精品美女| 亚洲天堂色综合久久| 久久热在线免费观看| 蜜乳视频一区二区三区| 91超碰国产在线观看| 无码人妻丰满熟妇区五路| 国内销魂老女人老泬| 一区二区三区资源视频| 色视频在线播放免费观看| 国产福利一区二区三区在线观看 | 在线免费观看a视频免费| 97香蕉久久国产超碰| 青青国产95免看视频| 伊人精品成人综合网| 超级黄肉动漫在线观看| av天堂a亚洲va天堂va里番| 亚洲在线免费观看18| 欧美三区四区在线视频| 伊人情人成综合视频| 日本东京热视频欧美视频| 国产91精品福利系列| 夜色17s精品人妻熟女av| 中文字幕福利视频在线一区| 欧美黑人性猛交小矮人| 熟妇高潮久久久久久久| 黄版视频在线免费观看| 成年人免费黄色av| 开心五月综合激情婷婷| 伊人网在线观看 视频一区| 亚洲国产美女主播在线观看| aaaa级少妇高潮在线观看| 91色老久久精品偷偷蜜臀| 韩国资源视频一区二区三区 | 91色老久久精品偷偷蜜臀| xxoo福利视频导航| 97成人老师在线视频| 国产午夜在线播放视频| 欧美vr专区日韩vr专区| 内地精品毛片在线观看| 国产福利小视频在线观看网站| 户外露出视频在线观看| 高清国产美女a一级毛片| 欧美精品999不卡| 久久久久九九九九九12| 亚洲蜜桃久久久久久| 欧美成人少妇人妻精品| 9999久久久久老熟妇二区| 欧美强奸视频在线观看| 偷拍熟女大胆免费视频| 一区二区三区国产在线成人av| 男人电影天堂在线观看| 国产青青青青草免费在线视频| 青娱乐这里只有精品| 亚州av嫩草av极品在线观看| 老司机免费视频福利0| 亚洲一区二区三区四区入口| 2020国产激情视频在线观看| 欧美亚洲愉拍一区二区三区| 一二三四区国产在线观看| 亚洲唯美激情综合四射| ysl蜜桃色7425| 亚洲综合天堂av网站在线观看| 得得爱在线视频观看| 日韩免费黄色片在线观看| 在线 制服 中文字幕 日韩| 二十四小时日本高清在线观看| 午夜久久久久欠久久久久| 鸡巴在里面福利视频在线观看| 亚洲在线观看中文字幕av| 亚洲欧美国产人成在线| 国产激情视频在线观看的| 1级黄色片在线观看| 欧美久久蜜臀蜜桃资源吧| 国产不卡免费在线观看| 亚洲第一页欧美第一页| 欧美大鸡吧男操女啊啊啊视频| 亚洲国产精品 久久久| 亚洲欧美另类丝袜另类自拍| 啪啪啪网站免费在线看| 福利一二三在线视频观看| 日韩人妻中文字幕区| 色视频免费观看网址| 综合激情网,激情五月| 麻豆白洁少妇在线播放| 熟妇人妻av无码中文字幕| 国长拍拍视频免费孕妇| 中国精品人妻一区二区| 91福利高清在线播放| 日本一本午夜在线播放| 亚洲欧美小说中文字幕| 手机看片1024精品国产| 熟妇人妻av无码中文字幕| 日本电影一级人妻在线播放四区| 成人午夜高清福利视频| 欧美一级aaaaaaa片| 婷婷色综合五月天视频| 大香蕉伊人97在线| 久久久人妻免费视频| 国产av高清二区三区| 国产av在线免费视频| 美国十次了亚洲天堂网国产| 大香焦一道本一区二区三区| 99久久久久久久久久久久久| 成人资源中文在线观看| 熟女人妻aⅴ一区二区三| 精品国产久久久久午夜精品av| 日本福利片在线播放| 五十岁熟女高潮喷水| 天天干夜夜操91视频网站| 欧美一区日韩二区三区四区| 一看就是假奶的av| 18福利视频在线观看| 久久精品四虎夜夜拍拍拍| 杜达雄啪啪毛片视频| 亚洲一区二区精品在线播放| 日本熟妇乱妇熟色视频| 看女人大BB群伦交| 欧美亚洲精品色图网站| 抽插小穴啊啊啊视频| 懂色av之国产精品| av一区二区三区蜜桃| 妈妈的朋友中字在线免费观看| 911美女片黄在线观看| 美女网站视频久久精品| 黑人黄色免费一级av| 夜夜操夜夜爱夜夜摸| 成年男女免费视频网站无毒| 国产一级一国产一级毛片| 中文字幕av特黄毛片| 欧美日韩一区二区三区成人影院| 91美女在线观看视频| 日本高清激情乱一区二区三区| 亚洲少妇色小说综合| 中文字幕日本一二三区| 天天在线播放日韩av| 天天插天天操天天射天天干| 亚洲综合首页综合在线观看| 91色老久久精品偷偷蜜臀| 奇米网首页神马久久| 亚洲熟女人妻自拍在线视频| 美国男的操女孩的小嫩逼| 92午夜免费福利视频www| 色屁屁一区二区三区在线观看| 538欧美在线观看一区二区三区| 成熟了的熟妇毛茸茸| 久久国产半精品99精品国产| 欧美日本国产一区二区| 中文字幕国产一区在线视频| 国产igao激情在线视频入口| 亚洲欧美国产一本综合首页| 日韩三级黄色大片在线观看| 九九热在线精品播放| 污网址在线观看视频| 放荡人妻极品少妇全集| 黄色av日韩在线观看| 国产欧美福利在线观看| 55夜色66夜色亚洲精品| 97精品人妻免费视频| 国产精品视频网站污污污 | xxoo福利视频导航| 中文字幕在线免费观看成人| 超peng视频在线免费播放97| 中文字幕在线免费观看成人| 一二区二区不卡视频| 欧美男男在线观看视频网站| aa福利影视在线观看| 免费绝清毛片a在线播放| 69精品互换人妻4p| 天天操天天搞天天操| 丝袜美女诱惑佐佐三上| xxoo福利视频导航| 少妇被粗大的猛进69视频| 狠狠干狠狠操免费视频| 2019年中文字幕在线播放视频| 在线免费观看欧美小视频| 中文字幕人妻一区色偷偷久久| 亚洲中文字幕在线视频观看二区 | 白白色在线免费视频发布视频| 国产成人在线观看hd| 911精产国品一二三产区区| 精品国产久久久久午夜精品av| 久久久久久久精品乱码| 中国精品人妻一区二区| 精品精品精品精品精品污污污污| 亚洲欧美日韩中文视频| 欧美成人少妇人妻精品| 欧美黄色一区二区三区视频| 免费中文三级在线观看| 国产精品性感美女视频| 午夜精品久久久久久久精品乱码 | 亚洲人妻系列在线视频| 99热在线只有的精品| 欧美亚洲愉拍一区二区三区| 精品国产av虐杀两警花| 国语对白性爱三级片免费看| 日本国产亚洲欧美色综合| 亚洲精品国品乱码久久久久| 亚洲人妻系列在线视频| 一区二区三区午夜福利在线| 国产探花自拍亚洲av| 91 精品视频在线看| 91九色91在线视频| 亚洲国产精品 久久久| 果冻麻豆一区二区三区| 夜夜躁av麻豆男| 国产人妻777人伦精品hd超碰| 亚洲少妇视频在线观看| 91九色pony蝌蚪| 东京热日韩av影片| 伊人网在线欧美日韩在线| 欧美vs亚洲vs日韩| 手机看电影一区二区三区| 黑人大巨屌操美女逼| 亚洲av中文无码网站| 亚洲一区亚洲二区成人福利| av男人站在线观看| 涩涩黄片在线免费观看| 一区二区九日韩美女| 四季av人妻一区二区三区| 亚洲gay视频在线观看| 女女抠逼白虎白丝袜| 自拍偷拍亚洲综合第一页| 国产高清在线观看av| 99久久久久久久久久久久久| 51vv精品视频在线观看| 91麻豆精品国产在线| 美国男的操女孩的小嫩逼| 在线观看网站伊人网| 干逼又爽又黄又免费的视频| 亚洲一区二区三区无码在线| 亚洲全国精品女人久久久| 福利一二三在线视频观看| 美女福利视频一区二区三区四区| 国产白丝一区二区三区av| 黄色av网址在线播放| 国产 亚洲 欧美 自拍| 一区二区三区观看在线| av人摸人人人澡人人超碰小说| 男女爱爱好爽视频免费看| 中文字幕av人妻一区二区三区| 亚洲无码专区中文字幕专区| 黄版视频在线免费观看| 中文字幕人妻精品精品| 日韩人妻一区二区三区在线观看| 激情九月天在线视频| 久久sm人妻中出精品一区二区| 亚洲天堂男人的天堂| 午夜3p福利视频合集| ass亚洲熟女ass| 欧美黄色性视频网站| 呻吟求饶的人妻中文字幕| 午夜亚洲国产精品中字| 西野翔人妻中文字幕中字在| 九九热视频1这里只有精品| 精品一区二区三区免费毛片W| 国产女人18毛片水真多精选| 亚洲国产电影的一区| 久久99国产中文丝袜| 亚洲一级熟妇丰满的女人| 老司机在线视频福利观看| 精品国产污污污污免费观看| 一区二区三区四区影片| 18福利视频在线观看| 女生抠逼自慰啊啊啊啊啊啊啊下载 | 欧美黑人性猛交小矮人| 中文字幕日韩人妻在线三区| av 资源在线播放| 91精产国品一二三产区区别网站| 大香蕉在线欧美在线视频| 国产经典精品欧美日韩| 亚洲欧美日韩中文在线观看| 午夜久久人妻一级内射av网址| 精产国品一二三77777| 青娱乐不卡视频在线| 丰满少妇_区二区三区| 最近最新最好看的中文字幕| 亚洲第一页欧美第一页| 豆豆专区操逼性视频在线| 久久99精品久久久久久三级| 国产探花自拍亚洲av| 亚洲美女a级黄色在线播放| 国产黄色主播网址大全在线播放| 欧美情色av在线观看| 亚洲字幕一区二区夜色av| 国产女人18毛片水真多精选| 国产午夜在线播放视频| 911美女片黄在线观看| 亚洲午夜熟女在线观看| 日韩人妻中文字幕区| 五月在线视频免费播放91| 大秀成年人国产精品视频 | 日本不卡视频一二三区| 一区二区三区观看在线| 精品美女洗澡一区二区| 精品不卡一区二区三区| 国产精美视频精品视频精品| 在线 制服 中文字幕 日韩| av激情四射五月婷婷| 欧美日本在线免费视频| 91福利高清在线播放| 日本东京热最新中文字幕| av激情四射五月婷婷| av网页免费在线观看| 欧美日韩一区二区三区成人影院| 99re这里是国产精品首页| 欧美在线观看一区二区不卡| 亚洲欧美不卡专业视频| 超碰在线免费观看视频97| 99久久99九九九99九| 可以直接看av网站| 日韩av水蜜桃一区二区三区| 国产,亚洲,欧美综合| 男人的天堂aⅴ在线| 国产欧美福利在线观看| 亚洲人人爽人人澡起碰av| 亚洲黑人欧美二区三区| 91精品一区一区三区| 最近中文字幕免费视频一| 欧美成人短视频在线播放| 91色乱一区二区三区| 69精品互换人妻4p| 91九色pony蝌蚪| 2020年亚洲男人天堂网| 69xx精品久久久久| 欧美成人性生活视频播放| 开心激情五月天作爱片| 日本清纯中文字幕版| 天天色天天射天天日天天干| 凹凸视频一区二区在线观看 | 欧美在线观看视频欧美| 豆豆专区操逼性视频在线| 亚洲午夜精品一级毛片app| 婷婷一区二区三区五月丁| 精品人妻人人做人人爽| 亚洲精品乱码久久久久app| 亚洲色大WWW永久网站| 黄色片黄色片黄色片黄色片黄色| 日本黄色一级电影网址| 色网站在线观看免费| 五十岁熟女高潮喷水| 久久久国产精品免费视频网| 不用付费特黄特色亚洲特级黄色片| 国产女主播在线观看一区| 亚洲宅男噜噜噜66在线观看| 夜夜躁av麻豆男| 亚洲欧美国产一本综合首页| 偷拍熟女大胆免费视频| 国内销魂老女人老泬| 亚洲成a人77777| 天天干夜夜撸天天操| 黄很色很在线免费视频网站| 亚洲男人天堂最新网址大全| 91精品国产欧美在线| 91精品麻豆91夜夜骚| 国产 亚洲 欧美 自拍| 性感人妻 中文字幕| 天天看天天爱天天日| 人妻视频网站快射视频网站| 老司机在线视频福利观看| 99久久国语露脸国产精品| 国产精品久久久99| 蜜乳av一区二区三区免费观看| 九一精品人妻一区二区三区| 国产精品成人免费电影| 综合激情网,激情五月| 51vv精品视频在线观看| 天天操天天干加勒比久久| 青青青在线视频观看97| 乌克兰美女操逼高清内射视频| 一区二区三区高清视频3| 大香蕉尹人在线最新| 91超碰九色porny| 天天操天天日天天插天天舔| 抽插小穴啊啊啊视频| 午夜精品老牛av一区二区三区| 丰满少妇人妻一区二区三区蜜桃| 538欧美在线观看一区二区三区 | 亚洲自拍偷拍一区二区中文字幕 | 网站在线观看蜜臀91| 狂操鸡巴小骚逼视频免费观看| 亚洲免费在线不卡视频| 全球高清中文字幕av| 福利视频免费在线播放| 不卡视频在线 欧美日韩| 黑鸡巴肏少妇逼视频| 最新日韩中文字幕免费在线观看| 国内销魂老女人老泬| 麻豆出品视频在线观看| 青娱乐免费视频一二三| 深夜福利免费观看在线看| 男人电影天堂在线观看| 国产一级一国产一级毛片 | 国产三级自拍视频在线观看网站 | 精品国产污污污污免费观看| 男人的天堂aⅴ在线| 欧美丝袜亚洲国产日韩| 国产男人的天堂一区| 天天日夜夜操人人爽| 最近最新最好看的中文字幕| 国产男女无套?免费网站下载| 亚洲国产日韩精品在线| 日本欧美国产在线一区| 亚洲乱码av一区二区蜜桃av| 亚洲综合另类欧美久久| 99久久国产精品免费消防器材| 大奶熟妇激情操逼逼| 亚洲a级视频在线播放| 68福利精品在线视频| 精久久久久久久久久久久| 欧美区日本区国产区| 国产精品网站亚洲发布| 91性高湖久久久久久久久久| 国产高清在线观看av| 丰满放荡熟妇在线播放| 自拍偷自拍亚洲精品10p| 亚洲制服丝袜网站中文字幕| 欧美老熟妇xxoo老妇| 天天干天天弄天天日| 色999日韩偷自拍拍免费| 国产白丝一区二区三区av| 男人av一区二区三区| 国产视频成人自拍蝌蚪视频| 欧美老熟妇xxoo老妇| 亚洲第一区av中文字幕| 69精品人妻久久久久久久久久久 | 男女69视频在线观看免费| 亚洲最强的25个城市| 日韩三级精品电影久久久久| 亚洲午夜高清在线观看| 日本老熟妇av老熟妇| 美利坚合众国av天堂| 手机看片福利一区二区三区四区| 久操资源在线免费播放| 久久99精品热在线观看| 无码精品黑人一区二区老人| 91青青青国产免费高清 | 午夜国产一区二区三区| 亚洲一区二区三区四区入口| 天天想要天天操天天干| 嗯~嗯~啊啊啊~高潮了软件| 亚洲欧美国产人成在线| 在线 激情 亚洲 视频| 大尺度久久久久久久| 亚洲国内精品久久久久久久| 黑人大吊大战亚洲女人。| 99热99这里免费的精品| 国产熟妇色xxⅹ交白浆视频| 神马不卡视频在线视频| 亚洲国产美女主播在线观看| www一区二区91| 天天夜夜久久精品综合| 先锋人妻啪啪中文字幕| 天天早上头和脸出汗是怎么办| 亚洲色视频在线播放网站| 五月婷婷伊人久久中文字幕| 91九色91在线视频| 欧美日本在线免费视频| 九色porny91国产| 男女69视频在线观看免费| 99福利一区二区视频| 日韩欧美一区二区三区免费看 | av里面的动作是真进去吗 | 大尺度av毛片在线网址| tushy一区二区三区视频| 欧美成人久久久桃色aa| 91精品国产91久久久久久密臀| 一区二区三区四区久久久久韩日| 欧美大胆a级视频秒播| 天天操天天舔天天做| 亚洲a区在线免费观看| 大鸡扒操大逼大片免费关看| 丝袜美女诱惑佐佐三上| 大陆中文字幕视频在线| 女女抠逼白虎白丝袜| 成人精品影视一区二区| 九九热精品视频在线播放| 九色porny91国产| 成人午夜av电影网| 成人精品动漫一区二区| 不卡高清一区二区三区| xxoo福利视频导航| 快使劲弄我视频在线播放| 亚洲成人av在线一区二区| 偷拍欧美日韩另类图片| 999久久久人妻精品一区| 九九热视频1这里只有精品| 日本熟妇乱妇熟色视频| 在线看日韩av不卡| 人人人妻人人人妻精品少妇| 老司机免费视频福利0| 亚洲黄色成人一级片| 五月的婷婷综合视频| 最新国产精品综合网高清| 老熟女xxxⅹhd老熟女性| 日本四十路人妻熟女| 国产中年夫妇激情高潮| 色哟哟亚洲乱码国产乱码精品精| 成人黄色录像在线观看| 最新福利二区三区视频| 中文字幕免费啪啪啪| av里面的动作是真进去吗| 1区3区4区产品乱入视频| 伊人精品久久一区二区 | 日本亚洲精品视频在线观看| 日本小视频一区二区| 99精品视频在线在线观看| 天天摸天天舔天天操天天日| 免费观看在线中文字幕视频| 欧美日韩综合精品无人区| 最新国产精品久久精品app| 中文字幕在线免费观看人妻| 美女av色播在线播放| 天天透天天舔天天操| 亚洲色视频在线播放网站| 日韩美精品成人一区二区三区四区| 大屁股熟女一区二区视频 | 91精品国产欧美在线| 久久久久夜色国产精品电影| 亚洲欧美激情国产综合久久久| 99久久国语露脸国产精品| 日本四十路人妻熟女| av天堂新资源在线| 亚洲色图日韩在线视频观看| 黄色大片一级老太太操逼| 91性高湖久久久久久久久久| 日韩无码国产一区二区| 久久中文字幕av一区二区| 日本熟妇乱妇熟色视频| 911精产国品一二三产区区| 老司机免费视频福利0| 亚洲美女色www色| 亚洲女人自熨在线视频| 国产成人在线观看视频播放| 亚洲男人的天堂最新网址| 69国产精品成人aaaaa片| 丰满少妇人妻一区二区三区蜜桃| 东京热日韩av影片| 国产熟女五十路一区二区三区| 抽插小穴啊啊啊视频| 亚洲国产精品一区二区第二页| 操死你美女在线视频| 七色福利视频在线观看| 精品国产无乱码一区二区三区| 亚洲自拍偷拍一区二区中文字幕 | 有码一区二区三区四区五区| 美国伦理片午夜理论片| 美女网站视频久久精品| 69精品互换人妻4p| 看女人大BB群伦交| 欧美精品999不卡| 亚洲欧美另类校园春色| 国产漂亮白嫩美女在线图片| 亚洲综合另类欧美久久| 美女露阴道让男人捅| 性感美女人妻久久久| 亚洲中文字幕无线乱码人妻精品 | 国产精美视频精品视频精品 | www国产亚洲精品久久久| 麻豆国产91制片厂| 青娱乐这里只有精品| 国产精品美女免费视频观看 | 国产一级一国产一级毛片| 天天干天天弄天天日| 天天在线播放日韩av| 国产做A爱免费视频在线观看| 九色91操最新在线观看网址| 4438x亚洲最大的成人| 日韩欧美一区二区三区免费看 | 亚洲 自拍 激情 另类| 中文字幕丰满子伦无码专区| 亚洲最强的25个城市| 第一福利视频在线观看| 国产精品国产三级在线高清观看| 最新国产精品拍在线观看| 欧美一级特黄大片做受99| 99re这里是国产精品首页| 青青操久久综合激情| 黄版视频在线免费观看| 日本人妻少妇xxxxxxx| 国产免费久久精品99re丫丫| 午夜精品久久秘?18免费观看| 伊人网在线观看 视频一区| 2020国产激情视频在线观看| 99久久人人爽亚洲精品美女 | 人妻色综合aaaaaa网| 91超碰九色porny| 日本少妇精品免费视频| 国产精品网站亚洲发布| 日韩黄色在线观看网站上| 极品少妇高潮喷水日出白浆| 亚洲美女露隐私av一区二区精品| 亚洲成人 国产精品| 18在线观看免费观看| av天堂新资源在线| 色999日韩偷自拍拍免费| 一区二区三区四区视频精品免费| 99999久久久精品| caopeng97在线观看视频| 强乱人妻中文字幕日本| 老鸭窝在线毛片观看免费播放| 在线视频自拍第三页| av无限看熟女人妻另类av| 青娱乐免费视频一二三| 青青青国产精品视频| 亚洲欧美另类校园春色| 日本少妇人妻中文在线| 久久久久九九九九九12| 伊人免费观看视频一| 天天爱天天日天天爽| 国产成人91色精品免费看片| 久久久久久久岛国免费观看| 亚洲国产美女主播在线观看| 欧美一区二区播放视频| 国产青青青青草免费在线视频| 一区二区三区av免费天天看| 51精品视频在线免费观看| 欧美极品少妇高潮喷水| 人妻色综合aaaaaa网| 制服丝袜 中文字幕 日韩| 亚洲色图日韩在线视频观看| avtt中文字幕手机版| 琪琪日本福利伦理视频| a级黄片免费观看| 欧美一区二区三区视频看| 亚洲美女a级黄色在线播放| 久草视频在线视频在线视频| 得得爱在线视频观看| 夜色17s精品人妻熟女av| 性色蜜桃臀x88av天美传媒| 亚洲gay视频在线观看| 每日更新日韩欧美在线| 亚洲AV无码一二三四区在线播放| 日日躁夜夜躁狠狠操| 97视频538在线观看| 欧美极品少妇高潮喷水| 国产探花自拍亚洲av| 亚洲高清一区二区三区久久| 欧美一级特黄大片在线| 免费中文字幕a级激情| 国长拍拍视频免费孕妇| 不用付费特黄特色亚洲特级黄色片| 最新激情中文字幕视频| 2021国产在线视频| 松本菜奈实最新av在线| 女人高潮潮呻吟喷水网站| 松本菜奈实最新av在线| 欧美成人红桃视频在线观看| 日韩av熟妇在线观看| 亚洲一区亚洲二区成人福利| 操烂你的骚逼天天欧美| 69av精品国产探花| 午夜免费福利老司机| 91佛爷视频在线观看| 男人电影天堂在线观看| 亚洲精品综合欧美精品综合| 久久久亚洲熟女一区二区| 自拍偷拍 国产激情| 精品人妻人人做人人爽| 99免费观看在线视频| 91大神在线免费观看视频| 上床啪啪啪免费视频| 黑人黄色免费一级av| 国产自拍偷拍视频在线免费观看 | 一区二区三区国产精华液区别大吗 | a级黄片免费观看| 美国伦理片午夜理论片| 五月天天堂视频在线| 99热在线只有的精品| 亚洲激情视频在线观看免费| 视频自拍偷拍视频自拍| 亚洲a区在线免费观看| 夜夜操夜夜爱夜夜摸| 亚洲同性同志一二三专区| 日韩三级精品电影久久久久 | 国产精品美女免费视频观看| 日韩在线 中文字幕| 久久久国产精品免费视频网| 中文字幕观看中文字幕免费 | av男人站在线观看| 日韩成人免费观看电影| 超peng视频在线免费播放97| 精品国产人伦一区二区三区| 日本熟妇乱妇熟色视频| 2020精品视频在线| 国产午夜在线播放视频| 欧美大胆a级视频秒播| 裸日本资源在线午夜| 美国男的操女孩的小嫩逼| 久久久久久免费观看av| 69视频在线精品国自产拍| 日本一区二区高清av中文| 人妻中文字幕亚洲在线| 日本一本午夜在线播放|