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

24小時熱門版塊排行榜    

查看: 9274  |  回復(fù): 57

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
引用回帖:
33樓: Originally posted by shenrenren at 2014-10-31 06:13:25

謝謝參與

[ 發(fā)自小木蟲客戶端 ]
不要說話
41樓2014-10-31 22:59:25
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
引用回帖:
36樓: Originally posted by magicmonk at 2014-10-31 08:32:55

謝謝參與

[ 發(fā)自小木蟲客戶端 ]
不要說話
42樓2014-10-31 22:59:44
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
引用回帖:
38樓: Originally posted by 子子丘水 at 2014-10-31 13:43:16

謝謝參與

[ 發(fā)自小木蟲客戶端 ]
不要說話
43樓2014-10-31 23:00:13
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
引用回帖:
40樓: Originally posted by gengxin60 at 2014-10-31 13:53:47

謝謝參與

[ 發(fā)自小木蟲客戶端 ]
不要說話
44樓2014-10-31 23:00:35
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
引用回帖:
39樓: Originally posted by EEB110 at 2014-10-31 13:50:49

謝謝參與

[ 發(fā)自小木蟲客戶端 ]
不要說話
45樓2014-10-31 23:00:59
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

hn19870519

木蟲 (正式寫手)
★ ★ ★ ★ ★ ★
小木蟲: 金幣+0.5, 給個紅包，謝謝回帖
ljw4010: 金幣+5, 謝謝交流！ 2014-11-11 09:50:01
ljw4010: 回帖置頂 2014-11-13 09:28:33
樓主開的這個帖是個好帖，不過被水軍拿走不少金幣～
我有個問題，在計(jì)算二維體系的時候，vasp計(jì)算的時候真空層如何處理。vasp計(jì)算的真空層和boltztrap的structr文件的真空層大小是否要保持一致？vasp計(jì)算的時候，如何測試K點(diǎn)對輸運(yùn)性質(zhì)的收斂性，我表示官方的手冊看不懂，寫的太簡單了。
最后一個問題：如果energy文件中減去了費(fèi)米能，那在intrance中是否還需要設(shè)置費(fèi)米能？
46樓2014-11-10 09:45:25
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
ljw4010: 回帖置頂 2014-11-15 21:43:43
引用回帖:
46樓: Originally posted by hn19870519 at 2014-11-10 09:45:25
樓主開的這個帖是個好帖，不過被水軍拿走不少金幣～
我有個問題，在計(jì)算二維體系的時候，vasp計(jì)算的時候真空層如何處理。vasp計(jì)算的真空層和boltztrap的structr文件的真空層大小是否要保持一致？vasp計(jì)算的時候，如 ...

你好，歡迎參與交流，我沒使用過VASP，也沒處理過二維體系！
1.幫你找了下關(guān)于真空層處理的一些問題,如下：
http://www.gaoyang168.com/html/201005/2032833.html
http://www.gaoyang168.com/html/201308/6236397.html
http://www.gaoyang168.com/html/201109/3601489.html
http://vasp1860.blog.edu.cn/home ... g&id=159735
2.收斂性問題：
http://www.gaoyang168.com/html/201311/6656831.html
http://www.gaoyang168.com/html/201212/5264940.html
3.使用Boltztrap時Energy文件有必要減去費(fèi)秘能么？如果減去了費(fèi)秘能，那只能把intrans中費(fèi)秘能設(shè)為0了
希望對你有幫助！
不要說話
47樓2014-11-11 09:49:12
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

hn19870519

木蟲 (正式寫手)
引用回帖:
47樓: Originally posted by ljw4010 at 2014-11-11 09:49:12
你好，歡迎參與交流，我沒使用過VASP，也沒處理過二維體系！
1.幫你找了下關(guān)于真空層處理的一些問題,如下：
http://www.gaoyang168.com/html/201005/2032833.html
http://www.gaoyang168.com/html/201308/6236397.html
http://em ...

非常感謝！
48樓2014-11-11 09:56:04
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

binshao1991

新蟲 (小有名氣)
★ ★ ★ ★ ★ ★
小木蟲: 金幣+0.5, 給個紅包，謝謝回帖
ljw4010: 金幣+5, 謝謝參與！ 2014-11-21 08:35:48
ljw4010: 回帖置頂 2014-11-21 08:58:40
樓主，能不能把給個VASP和boltztrap結(jié)合起來操作的例子��？很多人都是在用vasp算，本人是新鳥，極度期盼樓主回復(fù)
49樓2014-11-20 23:38:54
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖

ljw4010

榮譽(yù)版主 (職業(yè)作家)

小木蟲從頭派教主
ljw4010: 回帖置頂 2014-11-21 08:58:48
ljw4010: 取消置頂 2015-05-11 17:24:03
引用回帖:
49樓: Originally posted by binshao1991 at 2014-11-20 23:38:54
樓主，能不能把給個VASP和boltztrap結(jié)合起來操作的例子��？很多人都是在用vasp算，本人是新鳥，極度期盼樓主回復(fù)

你好，歡迎參與交流，由于本人主要是用WIEN2k,對vasp涉及少之又少，只能幫你找些參考資料，僅供參考！
1.VASP+BoltzTrap熱電計(jì)算程序安裝
1.準(zhǔn)備的文件：
（1）.VASP5.2.11程序包；
（2）.BoltzTrap_vasp 程序包（來源 Wenqing Zhang課題組）；

2.VASP編譯安裝
（1）.修改VASP的main.F
在 CALL DUMP_ALLOCATE(IO%IU6)前添加以下代碼，目的是產(chǎn)生SYSMETRY文件：
----------------------------------------------------------
!added by John. Yang for TE calculations
     OPEN(UNIT=38,FILE="SYMMETRY",STATUS="REPLACE"
     WRITE(38,"(I12)"NROTK
     DO NEDOS=1,NROTK   ! I DON'T WANT TO DEFINE A NEW INTEGER VARIABLE, SO HERE I USE NEDOS
        WRITE(38,"(3F10.5)" IGRPOP(1,1,NEDOS)*1.0,IGRPOP(1,2,NEDOS)*1.0,IGRPOP(1,3,NEDOS)*1.0
        WRITE(38,"(3F10.5)" IGRPOP(2,1,NEDOS)*1.0,IGRPOP(2,2,NEDOS)*1.0,IGRPOP(2,3,NEDOS)*1.0
        WRITE(38,"(3F10.5)" IGRPOP(3,1,NEDOS)*1.0,IGRPOP(3,2,NEDOS)*1.0,IGRPOP(3,3,NEDOS)*1.0
        WRITE(38,*)""
     ENDDO
     CLOSE(38)
-----------------------------------------------------------------------------
·        對vasp5.2.11版本，以上修改經(jīng)過測試，可編譯成功。
·        對vasp4.6版本，在"CALL TIMING"前添加，測試可編譯成功。

（2）編譯VASP,具體設(shè)置參照之前的文章：vasp5.2安裝(fftw3.3+Gotoblas)

3.BoltzTrap_vasp程序包編譯：
（1）參照組里的編譯環(huán)境，對Makefile進(jìn)行如下修改：主要改動BLAS,LAPACK,LIBS這幾個參數(shù)，以下僅顯示makefile的部分內(nèi)容；
------------------------------------------------------------
SHELL = /bin/sh
FC =ifort
FOPT = -FR -mp -C -g
FOPT = -FR -mp
LINKER  =       $(FC)
LFLAGS  =      
FGEN =
BLAS= /usr/local/lib/libgoto2.so
LAPACK= ../vasp5.2.11/vasp.5.lib/lapack_double.o
LIBS  = -L../vasp5.2.11/vasp.5.lib -ldmy \
    ../vasp5.2.11/vasp.5.lib/linpack_double.o $(LAPACK) \
    $(BLAS)
LDFLAGS = -L/export/mathlib/mkl/8.1.1/lib/em64t -Vaxlib -pthread
DESTDIR = .
EXECNAME = BoltzTrap_vasp
------------------------------------------------------------
（2）cp BoltzTrap_vasp ~/bin/

4,計(jì)算，
（1）VASP自洽計(jì)算（待考證是自洽計(jì)算還是DOS計(jì)算）
（2）建立文件：文件夾名.intrans
（摘自http://blog.sciencenet.cn/home.php?mod=space&uid=90975&do=blog&id=772892）
文件夾名.intrans(以下參數(shù)設(shè)置待考證，可以參考http://blog.sciencenet.cn/static/ueditor/dialogs/attachment/fileTypeImages/icon_default.pngUserGuide.pdf)
--------------------------------------------------------------
VASP                                                                        
0 0 0 0.0                          
0.233 0.0001 0.4  52
CALC                      # CALC (calculate expansion coeff), NOCALC read from file            
5                         # lpfac, number of latt-points per k                                 
BOLTZ                     # run mode (only BOLTZ is supported)                                 
0.15                       # (efcut) energy range of chemical potential                        
900. 900.                 # Tmax, temperature grid                                            
0.             # energyrange of bands given individual DOS output sig_xxx and dos_xxx (xxx is band number)
------------------------------------
（3）執(zhí)行：/home/zjma/setup/BoltzTrap_vasp/x_trans BoltzTrap_vasp
（4）結(jié)果在 文件夾名.trace及.condtens 文件中。

2.vasp計(jì)算的過程遇到的問題
（摘自http://blog.sina.com.cn/s/blog_b364ab230101e9dp.html）
最近在學(xué)vasp，這篇文章是百度文庫找到的，看了不錯，轉(zhuǎn)載一把。另外附上vasp程序，linux中下載后無須安裝即可使用。單機(jī)中可能會出現(xiàn)內(nèi)存溢出問題，可以放機(jī)群上使用。
01、第一原理計(jì)算的一些心得
（1）第一性原理其實(shí)是包括基于密度泛函的從頭算和基于Hartree-Fock自洽計(jì)算的從頭算，前者以電子密度作為基本變量（霍亨伯格-科洪定理），通過求解Kohn-Sham方程，迭代自洽得到體系的基態(tài)電子密度，然后求體系的基態(tài)性質(zhì)；后者則通過自洽求解Hartree-Fock方程，獲得體系的波函數(shù)，求基態(tài)性質(zhì)；
評述：K-S方程的計(jì)算水平達(dá)到了H-F水平，同時還考慮了電子間的交換關(guān)聯(lián)作用。
（2）關(guān)于DFT中密度泛函的Functional，其實(shí)是交換關(guān)聯(lián)泛函
包括LDA，GGA，雜化泛函等等
一般LDA為局域密度近似，在空間某點(diǎn)用均勻電子氣密度作為交換關(guān)聯(lián)泛函的唯一變量，多數(shù)為參數(shù)化的CA-PZ方案；
GGA為廣義梯度近似，不僅將電子密度作為交換關(guān)聯(lián)泛函的變量，也考慮了密度的梯度為變量，包括PBE,PW,RPBE等方案，BLYP泛函也屬于GGA；
此外還有一些雜化泛函，B3LYP等。
（3）關(guān)于贗勢
在處理計(jì)算體系中原子的電子態(tài)時，有兩種方法，一種是考慮所有電子，叫做全電子法，比如WIEN2K中的FLAPW方法(線性綴加平面波)；此外還有一種方法是只考慮價電子，而把芯電子和原子核構(gòu)成離子實(shí)放在一起考慮，即贗勢法，一般贗勢法是選取一個截?cái)喟霃剑財(cái)喟霃揭詢?nèi)，波函數(shù)變化較平滑，和真實(shí)的不同，截?cái)喟霃揭酝鈩t和真實(shí)情況相同，而且贗勢法得到的能量本征值和全電子法應(yīng)該相同。
贗勢包括模守恒和超軟，模守恒較硬，一般需要較大的截?cái)嗄埽泟輨t可以用較小的截?cái)嗄芗纯�。另外，模守恒勢的散射特性和全電子相同，因此一般紅外，拉曼等光譜的計(jì)算需要用模守恒勢。
    贗勢的測試標(biāo)準(zhǔn)應(yīng)是贗勢與全電子法計(jì)算結(jié)果的匹配度，而不是贗勢與實(shí)驗(yàn)結(jié)果的匹配度，因?yàn)楹蛯?shí)驗(yàn)結(jié)果的匹配可能是偶然的。
（4）關(guān)于收斂測試
（a）Ecut，也就是截?cái)嗄�，一般情況下，總能相對于不同Ecut做計(jì)算，當(dāng)Ecut增大時總能變化不明顯了即可；然而，在需要考慮體系應(yīng)力時，還需對應(yīng)力進(jìn)行收斂測試，而且應(yīng)力相對于Ecut的收斂要比總能更為苛刻，也就是某個截?cái)嗄芟驴偰芤呀?jīng)收斂了，但應(yīng)力未必收斂。
（b）K-point，即K網(wǎng)格，一般金屬需要較大的K網(wǎng)格，采用超晶胞時可以選用相對較小的K網(wǎng)格，但實(shí)際上還是要經(jīng)過測試。
（5）關(guān)于磁性
一般何時考慮自旋呢？舉例子，例如BaTiO3中，Ba、Ti和O分別為+2，+4和-2價，離子全部為各個軌道滿殼層的結(jié)構(gòu)，就不必考慮自旋了；對于BaMnO3中，由于Mn+3價時d軌道還有電子，但未滿，因此需考慮Mn的自旋，至于Ba和O則不必考慮。其實(shí)設(shè)定自旋就是給定一個原子磁矩的初始值，只在剛開始計(jì)算時作為初始值使用，具體的可參照磁性物理。
（6）關(guān)于幾何優(yōu)化
包括很多種了，比如晶格常數(shù)和原子位置同時優(yōu)化，只優(yōu)化原子位置，只優(yōu)化晶格常數(shù)，還有晶格常數(shù)和原子位置分開優(yōu)化等等。
在PRL一篇文章中見到過只優(yōu)化原子位置，晶格常數(shù)用實(shí)驗(yàn)值的例子（PRL 100, 186402 (2008)）；也見到過晶格常數(shù)先優(yōu)化，之后固定晶格常數(shù)優(yōu)化原子位置的情況；更多的情況則是Full geometry optimization。
一般情況下，也有不優(yōu)化幾何結(jié)構(gòu)直接計(jì)算電子結(jié)構(gòu)的，但是對于缺陷形成能的計(jì)算則往往要優(yōu)化。
（7）關(guān)于軟件
軟件大致分為基于平面波的軟件，如CASTEP、PWSCF和ABINIT等等，計(jì)算量大概和體系原子數(shù)目的三次方相關(guān)；還有基于原子軌道線性組合的軟件(LCAO)，比如openmx，siesta，dmol等，計(jì)算量和體系原子數(shù)目相關(guān)，一般可模擬較多原子數(shù)目的體系。
VASP是使用贗勢和平面波基組，進(jìn)行從頭量子力學(xué)分子動力學(xué)計(jì)算的軟件包，它基于CASTEP 1989版開發(fā)。VAMP/VASP中的方法基于有限溫度下的局域密度近似（用自由能作為變量）以及對每一MD步驟用有效矩陣對角方案和有效Pulay混合求解瞬時電子基態(tài)。這些技術(shù)可以避免原始的Car-Parrinello方法存在的一切問題，而后者是基于電子、離子運(yùn)動方程同時積分的方法。離子和電子的相互作用超緩Vanderbilt贗勢(US-PP)或投影擴(kuò)充波(PAW)方法描述。兩種技術(shù)都可以相當(dāng)程度地減少過渡金屬或第一行元素的每個原子所必需的平面波數(shù)量。力與張量可以用VAMP/VASP很容易地計(jì)算，用于把原子衰減到其瞬時基態(tài)中。
02、VASP程序的亮點(diǎn)：
1. VASP使用PAW方法或超軟贗勢，因此基組尺寸非常小，描述體材料一般需要每原子不超過100個平面波，大多數(shù)情況下甚至每原子50個平面波就能得到可靠結(jié)果。
2. 在平面波程序中，某些部分代碼的執(zhí)行是三次標(biāo)度。在VASP中，三次標(biāo)度部分的前因子足可忽略，導(dǎo)致關(guān)于體系尺寸的高效標(biāo)度。因此可以在實(shí)空間求解勢的非局域貢獻(xiàn)，并使正交化的次數(shù)最少。當(dāng)體系具有大約2000個電子能帶時，三次標(biāo)度部分與其它部分可比，因此VASP可用于直到4000個價電子的體系。
3. VASP使用傳統(tǒng)的自洽場循環(huán)計(jì)算電子基態(tài)。這一方案與數(shù)值方法組合會實(shí)現(xiàn)有效、穩(wěn)定、快速的Kohn-Sham方程自洽求解方案。程序使用的迭代矩陣對角化方案（RMM-DISS和分塊Davidson）可能是目前最快的方案。
4. VASP包含全功能的對稱性代碼，可以自動確定任意構(gòu)型的對稱性。
5. 對稱性代碼還用于設(shè)定Monkhorst-Pack特殊點(diǎn)，可以有效計(jì)算體材料和對稱的團(tuán)簇。Brillouin區(qū)的積分使用模糊方法或四面體方法。四面體方法可以用Blöchl校正去掉線性四面體方法的二次誤差，實(shí)現(xiàn)更快的k點(diǎn)收斂速度。
03、VASP 5.2的新功能：
1. 大規(guī)模并行計(jì)算需要較少的內(nèi)存。
2. 加入新的梯度校正泛函AM05和PBEsol；用標(biāo)準(zhǔn)PBE POTCAR文件提供新泛函；改善了單中心處理。
3. 離子位置和格矢中加入有限差分，從而得到二階導(dǎo)，用于計(jì)算原子間力常數(shù)和聲子（需要超晶胞近似），和彈性常數(shù)。計(jì)算中自動考慮對稱性。
4. 離子位置和靜電場中加入線性響應(yīng)，從而得到二階導(dǎo)，用于計(jì)算原子間力常數(shù)和聲子（需要超晶胞近似），Born有效電荷張量，靜態(tài)介電張量（電子和離子貢獻(xiàn)），內(nèi)應(yīng)變張量，壓電張量（電子和離子貢獻(xiàn)）。線性響應(yīng)只能用于局域和半局域泛函。
5. 精確的非局域交換和雜化泛函：Hartree-Fock方法；雜化泛函，特別是PBE0和HSE06；屏蔽交換；（實(shí)驗(yàn)性的）簡單模型勢GW-COHSEX，用于經(jīng)驗(yàn)的屏蔽交換內(nèi)核；（實(shí)驗(yàn)性的）雜化泛函B3LYP。
6. 通過本征態(tài)求和計(jì)算含頻介電張量：使用粒子無關(guān)近似，或通過GW的隨機(jī)相近似�？捎糜诰钟颍刖钟�，雜化泛函，屏蔽交換，和Hartree-Fock。
7. 完全含頻GW，速度達(dá)到等離子極點(diǎn)模型：單發(fā)G0W0；在G和W中迭代本征矢直至自洽；（實(shí)驗(yàn)性的）迭代G（也可以選W）本征矢的自洽GW；（實(shí)驗(yàn)性的）對相關(guān)能使用RPA近似的GW總能量；用LDA計(jì)算G和W的頂點(diǎn)校正（局域場效應(yīng)），僅能用于非自旋極化的情況；（實(shí)驗(yàn)性的）W的多體頂點(diǎn)校正，僅能用于非自旋極化的情況。
8. 實(shí)驗(yàn)性的功能：用TD-HF和TD-雜化泛函求解Cassida方程（僅能用于非自旋極化的Tamm-Dancoff近似）；GW頂點(diǎn)的Bethe-Salpeter（僅能用于非自旋極化的Tamm-Dancoff近似）。
1、VASP能夠進(jìn)行哪些過程的計(jì)算？怎樣設(shè)置？
我們平時最常用的研究方法是做單點(diǎn)能計(jì)算，結(jié)構(gòu)優(yōu)化、從頭計(jì)算的分子動力學(xué)和電子結(jié)構(gòu)相關(guān)性質(zhì)的計(jì)算。

一般我們的研究可以按照這樣的過程來進(jìn)行

如果要研究一個體系的最優(yōu)化構(gòu)型問題可以首先進(jìn)行結(jié)構(gòu)弛豫優(yōu)化，然后對優(yōu)化后的結(jié)構(gòu)進(jìn)行性質(zhì)計(jì)算或者單點(diǎn)能計(jì)算。

如果要研究一個體系的熱力學(xué)變化過程可以首先進(jìn)行分子動力學(xué)過程模擬，然后在某個溫度或壓強(qiáng)下進(jìn)行性質(zhì)計(jì)算或者單點(diǎn)能計(jì)算。

如果要研究一個體系的熱力學(xué)結(jié)構(gòu)變化可以首先在初始溫度下進(jìn)行NVT計(jì)算，然后進(jìn)行分子動力學(xué)退火，然后在結(jié)束溫度下進(jìn)行性質(zhì)計(jì)算研究。

2、什么是單點(diǎn)能計(jì)算(single point energy)？如何計(jì)算？

跟其它軟件類似，VASP具有單點(diǎn)能計(jì)算的功能。也就是說，對一個給定的固定不變的結(jié)構(gòu)（包括原子、分子、表面或體材料）能夠計(jì)算其總能，即靜態(tài)計(jì)算功能。

單點(diǎn)能計(jì)算需要的參數(shù)最少，最多只要在KPOINTS文件中設(shè)置一下合適的K點(diǎn)或者在INCAR文件中給定一個截?cái)嗄蹺NCUT就可以了。還有一個參數(shù)就是電子步的收斂標(biāo)準(zhǔn)的設(shè)置EDIFF，默認(rèn)值為EDIFF=1E-4，一般不需要修改這個值。

具體來說要計(jì)算單點(diǎn)能，只要在INCAR中設(shè)置IBRION=-1也就是讓離子不移動就可以了。

3、什么是結(jié)構(gòu)優(yōu)化(structure optimization)？如何計(jì)算？

結(jié)構(gòu)優(yōu)化又叫結(jié)構(gòu)弛豫（structure relax），是指通過對體系的坐標(biāo)進(jìn)行調(diào)整，使得其能量或內(nèi)力達(dá)到最小的過程，與動力學(xué)退火不同，它是一種在０Ｋ下用原子間靜力進(jìn)行優(yōu)化的方法�？梢哉J(rèn)為結(jié)構(gòu)優(yōu)化后的結(jié)構(gòu)是相對穩(wěn)定的基態(tài)結(jié)構(gòu)，能夠在實(shí)驗(yàn)之中獲得的幾率要大些（當(dāng)然這只是理論計(jì)算的結(jié)果，必須由實(shí)驗(yàn)來驗(yàn)證）。

一般要做弛豫計(jì)算，需要設(shè)置弛豫收斂標(biāo)準(zhǔn)，也就是告訴系統(tǒng)收斂達(dá)成的判據(jù)（convergence break condition)，當(dāng)系統(tǒng)檢測到能量變化減小到一個確定值時例如EDIFFG=1E-3時視為收斂中斷計(jì)算，移動離子位置嘗試進(jìn)行下一步計(jì)算。EDIFFG這個值可以為負(fù)，例如EDIFFG=-0.02，這時的收斂標(biāo)準(zhǔn)是當(dāng)系統(tǒng)發(fā)現(xiàn)所有離子間作用力都小于給定的數(shù)值，如０.０２eV/A時視為收斂而中斷。

弛豫計(jì)算主要有兩種方式：準(zhǔn)牛頓方法(quasi-Newton RMM-DIIS)和共軛梯度法(CG)兩種。準(zhǔn)牛頓方法計(jì)算速度較快，適合于初始結(jié)構(gòu)與平衡結(jié)構(gòu)（勢能面上全局最小值）比較接近的情況，而CG方法慢一些，找到全局最小的可能性也要大一些。選擇方法為IBRION=1時為準(zhǔn)牛頓方法而IBRION=2時為CG方法。

具體來說要做弛豫計(jì)算，設(shè)置IBRION=1或者2就可以了，其它參數(shù)根據(jù)需要來設(shè)置。NSW是進(jìn)行弛豫的最大步數(shù)，例如設(shè)置NSW=100，當(dāng)計(jì)算在100步之內(nèi)達(dá)到收斂時計(jì)算自動中斷，而100步內(nèi)沒有達(dá)到收斂的話系統(tǒng)將在第100步后強(qiáng)制中止（平常計(jì)算步數(shù)不會超過100步，超過100步可能是計(jì)算的體系出了問題）。參數(shù)通�？梢詮奈墨I(xiàn)中發(fā)現(xiàn)，例如收斂標(biāo)準(zhǔn)EDIFFG等。

有的時候我們需要一些帶限制條件的弛豫計(jì)算，例如凍結(jié)部分原子、限制自旋的計(jì)算等等。凍結(jié)部分原子可以在POSCAR文件中設(shè)置selective dynamic來實(shí)現(xiàn)。自旋多重度限制可以在INCAR中以NUPDOWN選項(xiàng)來設(shè)置。另外ISIF選項(xiàng)可以控制弛豫時的晶胞變化情況，例如晶胞的形狀和體積等。

費(fèi)米面附近能級電子分布的smearing是一種促進(jìn)收斂的有效方法，可能產(chǎn)生物理意義不明確的分?jǐn)?shù)占據(jù)態(tài)情況，不過問題不大。在INCAR文件中以ISMEAR來設(shè)置。一般來說K點(diǎn)只有一兩個的時候采用ISMEAR=0，金屬體材料用ISMEAR=1或2,半導(dǎo)體材料用ISMEAR=-5等等。不過有時電子步收斂速度依然很慢，還需要設(shè)置一些算法控制選項(xiàng)，例如設(shè)置ALGO=Very_Fast，減小真空層厚度，減少K點(diǎn)數(shù)目等。

弛豫是一種非常有效的分析計(jì)算手段，雖然是靜力學(xué)計(jì)算但是往往獲得一些動力學(xué)得不到的結(jié)果。

4、vasp的分子動力學(xué)模擬

vasp做分子動力學(xué)的好處，由于vasp是近些年開發(fā)的比較成熟的軟件，在做電子scf速度方面有較好的優(yōu)勢。缺點(diǎn)：可選系綜太少。

盡管如此，對于大多數(shù)有關(guān)分子動力學(xué)的任務(wù)還是可以勝任的。主要使用的系綜是NVT 和NVE。

一般做分子動力學(xué)的時候都需要較多原子，一般都超過100個。

當(dāng)原子數(shù)多的時候，k點(diǎn)實(shí)際就需要較少了。有的時候用一個k點(diǎn)就行，不過這都需要嚴(yán)格的測試。通常超過200個原子的時候，用一個k點(diǎn)，即Gamma點(diǎn)就可以了。

INCAR：
EDIFF   一般來說，用1E-4 或者1E-5都可以，這個參數(shù)只是對第一個離子步的自洽影響大一些，對于長時間的分子動力學(xué)的模擬，精度小一點(diǎn)也無所謂，但不能太小。
IBRION=0 分子動力學(xué)模擬
IALGO=48 一般用48，對于原子數(shù)較多，這個優(yōu)化方式較好。
NSW=1000   多少個時間步長。
POTIM=3  時間步長,單位fs, 通常1到3.
ISIF=2  計(jì)算外界的壓力.
NBLOCK= 1  多少個時間步長，寫一次CONTCAR，CHG和CHGCAR，PCDAT.
KBLOCK=50 NBLOCK*KBLOCK 個步長寫一次XDATCAR.（個離子步寫一次PCDAT.）
ISMEAR=-1  費(fèi)米迪拉克分布.
SIGMA =0.05 單位:電子伏
NELMIN=8  一般用6到8, 最小的電子scf數(shù).太少的話,收斂的不好.
LREAL=A
APACO=10 徑向分布函數(shù)距離, 單位是埃.
NPACO=200  徑向分布函數(shù)插的點(diǎn)數(shù).
LCHARG=F 盡量不寫電荷密度,否則CHG文件太大.
TEBEG=300  初始溫度.
TEEND=300 終態(tài)溫度。不設(shè)的話，等于TEBEG.
SMASS=-3  NVE ensemble;-1 用來做模擬退火。大于0 NVT 系綜。正確：SMASS=1,2,3 是沒有區(qū)別的。都是NVT ensemble。SMASS只要是大于0就是NVT系綜。

CONTCAR是每個離子步之后都會寫出來的，但是會用新的把老的覆蓋

CHG是在每10個離子步寫一次，不會覆蓋

CHGCAR是在任務(wù)正常結(jié)束之后才寫的。

5、收斂判據(jù)的選擇

結(jié)構(gòu)弛豫的判據(jù)一般有兩中選擇：能量和力。這兩者是相關(guān)的，理想情況下，能量收斂到基態(tài)，力也應(yīng)該是收斂到平衡態(tài)的。但是數(shù)值計(jì)算過程上的差異導(dǎo)致以二者為判據(jù)的收斂速度差異很大，力收斂速度絕大部分情況下都慢于能量收斂速度。這是因?yàn)榱Φ挠?jì)算是在能量的基礎(chǔ)上進(jìn)行的，能量對坐標(biāo)的一階導(dǎo)數(shù)得到力。計(jì)算量的增大和誤差的傳遞導(dǎo)致力收斂慢。

到底是以能量為收斂判據(jù)，還是以力為收斂判據(jù)呢？關(guān)心能量的人，覺得以能量為判據(jù)就夠了；關(guān)心力相關(guān)量的人，沒有選擇，只能用力作為收斂標(biāo)準(zhǔn)。對于超胞體系的結(jié)構(gòu)優(yōu)化，文獻(xiàn)大部分采用Gamma點(diǎn)做單點(diǎn)優(yōu)化。這個時候即使采用力為判據(jù)（EDIFFG=-0.02），在做靜態(tài)自洽計(jì)算能量的時候，會發(fā)現(xiàn)，原本已經(jīng)收斂得好好的力在不少敏感位置還是超過了結(jié)構(gòu)優(yōu)化時設(shè)置的標(biāo)準(zhǔn)。這個時候，是不是該懷疑對超胞僅做Gamma點(diǎn)結(jié)構(gòu)優(yōu)化的合理性呢？是不是要提高K點(diǎn)密度再做結(jié)構(gòu)優(yōu)化呢。

在我看來，這取決于所研究的問題的復(fù)雜程度。我們的計(jì)算從原胞開始，到超胞，到摻雜結(jié)構(gòu)，到吸附結(jié)構(gòu)，到反應(yīng)和解離。每一步都在增加復(fù)雜程度。結(jié)構(gòu)優(yōu)化終點(diǎn)與初始結(jié)構(gòu)是有關(guān)的，如果遇到對初始結(jié)構(gòu)敏感的優(yōu)化，那就頭疼了。而且，還要注意到，催化反應(yīng)不僅與原子本身及其化學(xué)環(huán)境有關(guān)，還會與幾何構(gòu)型有關(guān)。氣固催化反應(yīng)過程是電子的傳遞過程，也是分子拆分與重新組合的過程。如果優(yōu)化終點(diǎn)的構(gòu)型不同，可能會導(dǎo)致化學(xué)反應(yīng)的途徑上的差異。僅從這一點(diǎn)來看，第一性原理計(jì)算的復(fù)雜性，結(jié)果上的合理性判斷都不是手冊上寫的那么簡單。

對于涉及構(gòu)型敏感性的結(jié)構(gòu)優(yōu)化過程，我覺得，以力作為收斂判據(jù)更合適。而且需要在Gamma點(diǎn)優(yōu)化的基礎(chǔ)上再提高K點(diǎn)密度繼續(xù)優(yōu)化，直到靜態(tài)自洽計(jì)算時力達(dá)到收斂標(biāo)準(zhǔn)的。

6、結(jié)構(gòu)優(yōu)化參數(shù)設(shè)置

結(jié)構(gòu)優(yōu)化，或者叫弛豫，是后續(xù)計(jì)算的基礎(chǔ)。其收斂性受兩個主要因素影響：初始結(jié)構(gòu)的合理性和弛豫參數(shù)的設(shè)置

初始結(jié)構(gòu)

初始結(jié)構(gòu)包括原子堆積方式，和自旋、磁性、電荷、偶極等具有明確物理意義的模型相關(guān)參數(shù)。比如摻雜，表面吸附，空位等結(jié)構(gòu)，初始原子的距離，角度等的設(shè)置需要有一定的經(jīng)驗(yàn)積累。DFT計(jì)算短程強(qiáng)相互作用（相對于范德華力），如果初始距離設(shè)置過遠(yuǎn)（如超過4埃），則明顯導(dǎo)致收斂很慢甚至得到不合理的結(jié)果。

比較好的設(shè)置方法可以參照鍵長。比如CO在O頂位的吸附，可以參照CO2中C-O鍵長來設(shè)置（如增長20%）。也可以參照文獻(xiàn)。記住一些常見鍵長，典型晶體中原子間距離等參數(shù)，有助于提高初始結(jié)構(gòu)設(shè)置的合理性。實(shí)在不行，可以先在小體系上測試，然后再放到大體系中算。

弛豫參數(shù)

弛豫參數(shù)對收斂速度影響很大，這一點(diǎn)在計(jì)算工作沒有全部鋪開時可能不會覺察到有什么不妥，反正就給NSW設(shè)置個“無窮大”的數(shù)，最后總會有結(jié)果的。但是，時間是寶貴的，恰當(dāng)?shù)脑O(shè)置3小時就收斂的結(jié)果，不恰當(dāng)?shù)脑O(shè)置可能要一個白天加一個黑夜。如果你趕文章或者趕著畢業(yè)，你就知道這意味這什么。

結(jié)構(gòu)優(yōu)化分電子迭代和離子弛豫兩個嵌套的過程。電子迭代自洽的速度，有四個響很大的因素：初始結(jié)構(gòu)的合理性，k點(diǎn)密度，是否考慮自旋和高斯展寬（SIGMA）；離子弛豫的收斂速度，有三個很大的影響因素：弛豫方法（IBRION）,步長（POTIM）和收斂判據(jù)（EDIFFG）.

一般來說，針對理論催化的計(jì)算，初始結(jié)構(gòu)都是不太合理的。因此一開始采用很粗糙的優(yōu)化（EDIFF=0.001，EDIFFG=-0.2），很低的K點(diǎn)密度(Gamma)，不考慮自旋就可以了，這樣NSW<60的設(shè)置就比較好。其它參數(shù)可以默認(rèn)。

經(jīng)過第一輪優(yōu)化，就可以進(jìn)入下一步細(xì)致的優(yōu)化了。就我的經(jīng)驗(yàn)，EDIFF=1E-4,EDIFFG=-0.05，不考慮自旋，IBRION=2，其它默認(rèn)，NSW=100;跑完后可以設(shè)置IBRION = 1，減小OPTIM（默認(rèn)為0.5，可以設(shè)置0.2）繼續(xù)優(yōu)化。

優(yōu)化的時候讓它自己悶頭跑是不對的，經(jīng)�？纯粗虚g過程，根據(jù)情況調(diào)節(jié)優(yōu)化參數(shù)是可以很好的提高優(yōu)化速度。這個時候，提交兩個以上的任務(wù)排隊(duì)是好的方式，一個在調(diào)整的時候，下一個可以接著運(yùn)行，不會因?yàn)橥Ｏ庐?dāng)前任務(wù)導(dǎo)致機(jī)器空閑。

無論結(jié)構(gòu)優(yōu)化還是靜態(tài)自洽，電子步的收斂也常常讓新手頭痛。如果電子步不能在40步內(nèi)收斂，要么是參數(shù)設(shè)置的問題，要么是初始模型太糟糕（糟糕的不是一點(diǎn)點(diǎn)）。

靜態(tài)自洽過程電子步不收斂一般是參數(shù)設(shè)置有問題。這個時候，改變迭代算法（ALGO），提高高斯展寬（SIGMA增加）,設(shè)置自洽延遲（NELMDL）都是不錯的方法。對于大體系比較難收斂的話，可以先調(diào)節(jié)AMIN,BMIX跑十多步，得到電荷密度和波函數(shù)，再重新計(jì)算。實(shí)在沒辦法了，可以先放任它跑40步，沒有收斂的跡象的話，停下來，得到電荷密度和波函數(shù)后重新計(jì)算。一般都能在40步內(nèi)收斂。

對于離子弛豫過程，不調(diào)節(jié)關(guān)系也不大。開始兩個離子步可能要跑滿60步（默認(rèn)的），后面就會越來越快了。

總的說來，一般入門者，多看手冊，多想多理解，多上機(jī)實(shí)踐總結(jié)，比較容易提高到一個熟練操作工的水平。

如果要想做到“精確打擊”，做到能在問題始發(fā)的時候就立刻采取有效措施來解決，就需要回歸基礎(chǔ)理論和計(jì)算方法上來了。

7、優(yōu)化結(jié)果對初始結(jié)構(gòu)和“優(yōu)化路徑”的依賴

原子吸附問題不大，但是小分子吸附，存在初始構(gòu)型上的差異。slab上水平放置，還是垂直放置，可能導(dǎo)致收斂結(jié)果上的差異。根據(jù)H-K理論，理想情況下，優(yōu)化得到的應(yīng)該是全局最小，但在數(shù)值計(jì)算的時候可能經(jīng)常碰到不是全局最小的情況。實(shí)際操作中發(fā)現(xiàn)，多個不同初始結(jié)構(gòu)優(yōu)化收斂后在能量和結(jié)構(gòu)上存在一定差異。

為了加快收斂速度，特別是對于表面-分子吸附結(jié)構(gòu)，初始放松約束，比如EDIFF=1E-3，EDIFFG=-0.3,NSW=30可能是很好的設(shè)置。但是下面的情況應(yīng)當(dāng)慎重：

EDIFF=1E-3;

EDIFFG=-0.1;！或者更小

NSW=500；！或者更大

電子步收斂約束較小，而離子步約束偏大，離子步數(shù)又很多，這種情況下，可能導(dǎo)致的結(jié)果是結(jié)構(gòu)弛豫到嚴(yán)重未知的區(qū)間。

再在這個基礎(chǔ)上提高約束來優(yōu)化，可能就是徒勞的了——結(jié)果不可逆轉(zhuǎn)的偏向不正常的區(qū)間。

好的做法，是對初始結(jié)構(gòu)做比較松弛的約束，弛豫離子步NSW應(yīng)該限制在一個較小的數(shù)值內(nèi)。EDIFF=1E-3的話，EDIFFG也最好是偏大一些，如-0.3而不是-0.1. 這樣可以在較少的步數(shù)內(nèi)達(dá)到初步收斂。

對于遠(yuǎn)離基態(tài)的初始結(jié)構(gòu)，一開始在非常松弛的約束下跑若干離子步，時間上帶來的好處是很大的。對于100個原子的體系用vasp做Gamma點(diǎn)優(yōu)化，如果一開始就是正常優(yōu)化（EDIFF=1E-4,EDIFFG=-0.02）設(shè)置，開始十個離子步可能都要花上幾個小時。如果這個時候才發(fā)現(xiàn)輸入文件有錯誤，那下午的時間就白費(fèi)了，順便帶上晚上機(jī)器空轉(zhuǎn)。

所以，我習(xí)慣的做法，是在初始幾步優(yōu)化后，會用xcrysden 檢查一下XDATCAR中的數(shù)據(jù)，用xdat2xyz.pl生成movie.xyz，然后看看弛豫過程是不是按照設(shè)想的那樣。后續(xù)過程跑完一個收斂過程，就再檢查一下movie.xyz。如此這般，才放心的展開后續(xù)計(jì)算。

8、目的導(dǎo)向的結(jié)構(gòu)優(yōu)化

結(jié)構(gòu)優(yōu)化到這個階段，是高級的了。為了得到特定結(jié)構(gòu)，或者為了驗(yàn)證某些猜想，需要設(shè)計(jì)合理的初始結(jié)構(gòu)，然后在這個基礎(chǔ)上小心優(yōu)化，比如POTIM=0.1跑幾步看看，然后修改優(yōu)化參數(shù)。

我遇到過的一件跟結(jié)構(gòu)優(yōu)化關(guān)系很大的算例是CeO2氧空位結(jié)構(gòu)電子局域的問題(http://www.gaoyang168.com/bbs/viewthread.php?tid=2954558)。按照一般方式（從優(yōu)化好的bulk建slab模型，然后優(yōu)化）得到一個O空位留下的兩個電子均勻局域到O次外層三個Ce原子上，得到空位形成能2.34eV.經(jīng)高人指點(diǎn)后，調(diào)節(jié)空位附近O原子位置，打破對稱性后重新優(yōu)化，兩個電子完美的局域到兩個Ce原子上了。并且空位形成能降低到2.0X eV。從這個例子可以看到，結(jié)構(gòu)優(yōu)化存在不少技巧的，這些技巧建立在研究者對模擬對象的物理意義的理解上。對物理圖像的直觀深入理解，才能做好模型預(yù)設(shè)，在此引導(dǎo)下才可能有目的的優(yōu)化出不比尋常的結(jié)果。




目前第一性原理理論中的交換關(guān)聯(lián)泛函部分包含經(jīng)驗(yàn)參數(shù)�？紤]這一點(diǎn)對優(yōu)化結(jié)果的影響也很有意思。比如有專家提到，DFT+U參數(shù)對某些結(jié)構(gòu)的收斂終態(tài)構(gòu)型有影響。構(gòu)型的變化可能影響表面反應(yīng)過程�；谶@一點(diǎn)，一個好的計(jì)算研究可能就出來了。




真實(shí)過程總是復(fù)雜多變的。無論何種模擬，估計(jì)都可以找到一些試驗(yàn)現(xiàn)象來驗(yàn)證。但是到底應(yīng)該如何評判模擬結(jié)果，如何從第一性原理研究中得出有意義的結(jié)論需要很好的洞察力。這樣的模擬不見得就必須建立的試驗(yàn)的基礎(chǔ)上，完全憑空設(shè)計(jì)的模型有可能更能優(yōu)美的解釋本質(zhì)。




9、在單機(jī)上計(jì)算態(tài)密度好像不會出問題。我先談一下我的看法：
第一個WARNING，可以在INCAR文件中設(shè)置NGX,NGY和NGZ的值，設(shè)置的值要足夠大，就可以消除這個warning。設(shè)置多大合適呢？這就要用到編譯vasp時，同時也編譯得到的make param小程序，make param可以幫助你預(yù)先檢查你設(shè)置的文件是否正確，以及某些參數(shù)的值是否合適。要得到合適的NGX,NGY,NGZ以及NBANDS,先在INCAR中不設(shè)置這些參數(shù)的值，然后運(yùn)行makeparam >param.inc，其中param.inc是包含了輸出結(jié)果的文件，在param.inc文件中你可以看到這些參數(shù)的值，以及計(jì)算大概需要多少的內(nèi)存。然后把param.inc文件中的NGX,NGY,NGZ和NBANDS的值拷貝到INCAR文件中。
第二個是計(jì)算態(tài)密度時，我個人的做法是，一般把KPOINTS文件中的k點(diǎn)增多，然后把INCAR文件中的ISTART=1,ICHARG=11，當(dāng)然還設(shè)置RWIGS。最后把靜止自洽計(jì)算得到的CHG和CHGCAR文件拷貝到當(dāng)前目錄下。從我在單機(jī)上的計(jì)算來看，沒有WAVECAR文件也是可以計(jì)算態(tài)密度的。我想你出現(xiàn)的這個問題，可能是你cluster上計(jì)算時，每個節(jié)點(diǎn)上的CHGCAR和WAVECAR文件不一致造成的。
第三個是當(dāng)k點(diǎn)數(shù)增加了，會出現(xiàn)一個WARING，要把此WARNING消失掉，在INCAR文件中設(shè)置NELMDL，它的值小于等于默認(rèn)值（默認(rèn)值好像是-5,你可以設(shè)為-6)。沒有cluster的系統(tǒng)用來計(jì)算，也沒有這樣的經(jīng)歷，我僅從在單機(jī)上的計(jì)算經(jīng)驗(yàn)來談，有錯還請包涵
10、如何用VASP計(jì)算鐵磁、反鐵磁和順磁

順磁，意味進(jìn)行non-spin polarized的計(jì)算，也就是ISPIN=1。

鐵磁，意味進(jìn)行spin-polarized的計(jì)算，ISPIN=2，而且每個磁性原子的初始磁矩設(shè)置為一樣的值，也就是磁性原子的MAGMOM設(shè)置為一樣的值。對非磁性原子也可以設(shè)置成一樣的非零值（與磁性原子的一樣）或零，最后收斂的結(jié)果，非磁性原子的local磁矩很小，快接近0，很小的情況，很可能意味著真的是非磁性原子也會被極化而出現(xiàn)很小的local磁矩。

反鐵磁，也意味著要進(jìn)行spin-polarized的計(jì)算，ISPIN=2，這是需采用反鐵磁的磁胞來進(jìn)行計(jì)算，意味著此時計(jì)算所采用的晶胞不再是鐵磁計(jì)算時的最小原胞。比如對鐵晶體的鐵磁狀態(tài)，你可以采用bcc的原胞來計(jì)算，但是在進(jìn)行反鐵磁的Fe計(jì)算，這是你需要采用sc的結(jié)構(gòu)來計(jì)算，計(jì)算的晶胞中包括兩個原子，你要設(shè)置一個原子的MAGMOM為正的，另一個原子的MAGMOM設(shè)置為負(fù)，但是它們的絕對值一樣。因此在進(jìn)行反鐵磁的計(jì)算時，應(yīng)該確定好反鐵磁的磁胞，以及磁序，要判斷哪種磁序和磁胞是最可能的反鐵磁狀態(tài)，那只能是先做好各種可能的排列組合，然后分別計(jì)算這些可能組合的情況，最后比較它們的總能，總能最低的就是可能的磁序。同樣也可以與它們同鐵磁或順磁的進(jìn)行比較。了解到該材料究竟是鐵磁的、還是順磁或反鐵磁的。

亞鐵磁，也意味要進(jìn)行spin-polarized的計(jì)算，ISPIN=2，與反鐵磁的計(jì)算類似，不同的是原子正負(fù)磁矩的絕對值不是樣大。非共線的磁性，那需采用專門的non-collinear的來進(jìn)行計(jì)算，除了要設(shè)置ISPIN，MAGMOM的設(shè)置還需要指定每個原子在x,y,z方向上的大小。這種情況會復(fù)雜一些。

舉個例子來說，對于Mn-Cu(001)c(2x2)這種體系，原胞里面有2個Mn原子，那么你直接讓兩個Mn原子的MAGMOM的絕對值一樣，符號相反就可以了，再加上ISPIN=2。這樣就可以實(shí)現(xiàn)進(jìn)行反鐵磁的計(jì)算了

11、vasp在計(jì)算磁性的時候，oszicar中得到的磁矩和outcar中得到各原子磁矩之和不一致在投稿的是否曾碰到有審稿人質(zhì)疑，對于這個不一致你們一般是怎么解釋的了？

答：OSZICAR中得到的磁矩是OUTCAR中最后一步得到的總磁矩是相等的。總磁矩和各原子的磁矩(RMT球內(nèi)的磁矩)之和之差就是間隙區(qū)的磁矩。因?yàn)橛虚g隙區(qū)存在,不一致是正常的。

12、磁性計(jì)算應(yīng)該比較負(fù)責(zé)。你應(yīng)該還使用別的程序計(jì)算過磁性，與vasp結(jié)果比較是否一致，對磁性計(jì)算采用的程序有什么推薦。
ps：由于曾使用vasp和dmol算過非周期體系磁性，結(jié)構(gòu)對磁性影響非常大，因此使用這兩個程序計(jì)算的磁性要一致很麻煩。還不敢確定到底是哪個程序可能不可靠。
答：如果算磁性，全電子的結(jié)果更精確，我的一些計(jì)算結(jié)果顯示磁性原子對在最近鄰的位置時，PAW與FPLAW給出的能量差不一致，在長程時符合的很好。雖然并沒有改變定性結(jié)論。感覺PAW似乎不能很好地描述較強(qiáng)耦合。我試圖在找出原因，主要使用exciting和vasp做比較。計(jì)算磁性推薦使用FP-LAPW, FP-LMTO, FPLO很吸引人(不過是商業(yè)的)，后者是O(N)算法。
13、vasp 學(xué)習(xí)筆記POTCAR 的建立
POTCAR將要告訴vasp計(jì)算的系統(tǒng)中所包含的各種元素的贗勢pesudopotential，vasp本身就帶有比較完善的贗勢包，我們需要做的就是選擇我們需要具體哪種贗勢，然后把相應(yīng)的文件拷貝形成我們具體的POTCAR文件。我們以GaAs為例。
1）贗勢的選擇：
vasp的贗勢文件放在目錄~/vasp/potentials 下，可以看到該目錄又包含五個子目錄pot pot_GGA potpaw potpaw_GGA potpaw_PBE ，其中每一個子目錄對應(yīng)一種贗勢形式。
贗勢按產(chǎn)生方法可以分為PP (standard pesudopotential，其中大部分是USPP, ultrasoft pesudopotential) 和PAW (projector augmented wave method)。按交換關(guān)聯(lián)函數(shù)的不同又可以有LDA (local density approximation) 和GGA (generalized gradient approximation)，其中GGA之下又可以再分為PW91和PBE。
以上各個目錄對應(yīng)起來分別是pot ==> PP, LDA ; pot_GGA ==> PP, GGA ; potpaw ==> PAW, LDA ; potpaw_GGA ==> PAW, GGA, PW91 ; potpaw_PBE ==> PAW , GGA, PBE。選擇某個目錄進(jìn)去，我們還會發(fā)現(xiàn)對應(yīng)每種元素往往還會有多種贗勢存在。這是因?yàn)楦鶕?jù)對截?cái)嗄芰康倪x取不同還可以分為Ga,Ga_s,Ga_h，或者根據(jù)半芯態(tài)的不同還可以分為Ga,Ga_sv,Ga_pv的不同。
一般推薦選取PAW_PBE。其中各個元素具體推薦哪種形式的贗勢可以參考vasp workshop中有關(guān)贗勢部分的ppt。當(dāng)然自己能測試之后在選擇是最好不過的了，以后再聊。
2）.POTCAR的建立：
選好哪一種贗勢之后，進(jìn)入對應(yīng)的目錄，你會看到里邊有這么幾個文件，POTCAR.Z PSCTR.Z V_RHFIN.Z WS_FTP.LOG 。我們需要的是第一個。把它解壓，如zcat POTCAR.Z > Ga 。對As元素我們也可以類似得到一個As文件。用cp 命令或者mv 命令把這兩個文件都移到我們的工作目錄里。然后再用cat 命令把這兩個文件合并在一起，如cat Ga As > POTCAR ，這樣就得到了我們需要的POTCAR。同理，有多個元素的POTCAR也可以這樣產(chǎn)生。這里需要注意的是，記住元素的排列順序，以后在POSCAR里各個元素的排列就是按著這里來的。
3）.POTCAR里的信息：
如果你想看POTCAR長什么樣，可以用vim POTCAR 命令，進(jìn)去后可以用上下鍵移動光標(biāo)。想出來的時候，可以敲入:q!就可以。具體的vim的命令可以在網(wǎng)上查到。一般我會看POTCAR里的截?cái)嗄芰繛槎啻�，用grep -in "enmax" POTCAR 。
據(jù)說B3LYP的贗勢計(jì)算比較準(zhǔn)，我在MS上面測試過，好像DOS和能帶圖的計(jì)算確實(shí)比較準(zhǔn)。不過不知道vasp有沒有類似的贗勢包。
hybrid functional的計(jì)算，并不需要特定的hybrid functional 的贗勢。大部分就是基于GGA-PBE的贗勢來做，也就是芯電子與價電子的交換關(guān)聯(lián)作用，以及芯電子與芯電子的交換關(guān)聯(lián)作用還是基于GGA-PBE的，只是將價電子與價電子的交換關(guān)聯(lián)作用通過hybrid functional交換關(guān)聯(lián)來描述。
謝謝老師的解答。那具體操作是不是像網(wǎng)上寫的那樣，使用GGA的贗勢，設(shè)置GGA = B3，然后更改POTCAR里面的LEXCH =B3就行了。我試過了，可以跑，不過結(jié)果沒做詳細(xì)的分析。
14、VASP中所有能量的物理意義及它們之間的區(qū)別，讓你徹底搞清楚VASP的所有能量
（一）首先我們應(yīng)明白，固體的結(jié)合能就是固體的內(nèi)能E（結(jié)合）=U（內(nèi)能），
原因如下：
一般情況都把孤立原子的能量作為能量參考點(diǎn)。前段時間有個同學(xué)問VASP中得出的絕對能量是相對于什么的，其實(shí)就是相對孤立原子得。
（二）其次我們根據(jù)自由能與內(nèi)能之間的關(guān)系F=U-TS
而且我們都知道VASP的所有計(jì)算都是在絕對0度下的情況，T=0代入上式，有F=U。所以結(jié)合就等于內(nèi)能等于自由能�？隙ㄓ蠪ree energy TOTEN=energy without entropy恒成立...
這時候肯定有人會說不對啊，可以看VASP手冊，候博的參考書作證，肯定不對得。
現(xiàn)在我告訴你確實(shí)它們二者確實(shí)有區(qū)別，區(qū)別在下面的情況
（1）當(dāng)我們用ISMEAR=-5時，費(fèi)米能這兒沒有展寬，它算出來的就是完全在絕對0度的能量。Free energy TOTEN=energy without entropy恒成立。
（2）有時為了在數(shù)學(xué)上處理的方便，為了更容易積分，我們也用ISMEAR！=-5（！=是不等于的意思）的方法，這個時候費(fèi)米能這兒有一定的展寬。此時，我們?nèi)菀紫氲�，有展寬不就是相�?dāng)有一定的熵值嗎？所以這個時候雖然算的是絕對0度的情況，但是有一定的熵值（我們應(yīng)明白，這個熵值不是由一定的溫度帶來的，而是數(shù)學(xué)處理的結(jié)果）。所以在SMEAR！=-5的方法我們會發(fā)現(xiàn)Free energy TOTEN和energy without entropy有一定的差別。此時energy without entropy是Free energy TOTEN在SIGMA趨于0的極限。
注意:（1）有人在算單個原子的能量時會發(fā)現(xiàn)單個原子的能量雖然很小但并不是0,但是按我上面的推導(dǎo),固體中的結(jié)合能是相對孤立體系的能量而來的,所以單個原子得到的TOTEN肯定是0啊,原因在于我們的POTCAR不可能絕對合理,而且我們也知道計(jì)算單個原子的能量就是為了檢測贗勢,單原子得到的TOTEN越小說明贗勢越好。但一般不會正好是0.對這個說法我還存在點(diǎn)疑問，寫在了最后面。
（2）如果你注意的話，energy without entropy與Free energy TOTEN在SIGMA趨于0也不是完全相等，但是也會發(fā)現(xiàn)它們之間的差別在10E-3左右，原因在于計(jì)算機(jī)求積分、求極限不能像我們?nèi)艘粯舆_(dá)到任意的精度。
15、VASP中過渡態(tài)計(jì)算設(shè)置的一點(diǎn)體會
計(jì)算過渡態(tài)先要擺正心態(tài)，不急于下手。步驟如下：
（1）做模型，初態(tài)IS和終態(tài)FS，分別結(jié)構(gòu)優(yōu)化到基態(tài)；
（2）線形插入images: nebmake.pl POSCAR.IS POSCAR.FS N
N為image個數(shù)。
（3）nebmovie.pl,生成movie.xyz。用Xcrysden --xyz movie.xyz 反復(fù)觀看動畫，仔細(xì)檢查過程的合理性。這里要提醒，POSCAR.IS 和POSCAR.FS中原子坐標(biāo)列表的順序必須對應(yīng)。
（4）寫INCAR,選IOPT。注意，最好忘記vasp自帶的NEB，而全部改用包含vtstool的vasp. IBRION=3,POTIM=0關(guān)閉vasp自帶的NEB功能。
（5）過渡態(tài)計(jì)算第一個離子步最耗時，也最容易出問題，也是模型設(shè)計(jì)合理性檢驗(yàn)的首要環(huán)節(jié)。所以可以選小一些的ENCUT，可以不用考慮自旋(ISPIN=1)，也不用考慮DFT+U。而且用最快最粗糙的算法（IOPT=3,其他默認(rèn)）。
（6）帶vtstool的vasp-ClNEB(NEB)過渡態(tài)計(jì)算ICHAIN=0作為入口，這個也是默認(rèn)的。LCLIMB=TRUE也是默認(rèn)的。如果不要climb image,可以設(shè)置LCLIMB = False.
（7）收斂判據(jù)EDIFFG<0。過渡態(tài)計(jì)算要以力為收斂判據(jù)，而不是能量。一般EDIFFG=-0.05就可以接受，-0.02或者-0.01更好。但是作為開始的過渡態(tài)計(jì)算，可以設(shè)置很寬的收斂條件，如EDIFFG=-1.
（8）初步過渡態(tài)收斂后,修改INCAR中的優(yōu)化器（IOPT），并修改相應(yīng)參數(shù)（參考vtstool官方論壇），EDIFFG改�。ㄈ�-0.05），然后運(yùn)行vfin.pl，這個腳本自動幫你準(zhǔn)備在原來的基礎(chǔ)上繼續(xù)運(yùn)行新的過渡態(tài)計(jì)算（完成cp CONTCAR POSCAR, 保留電荷密度和波函數(shù)的操作）。
（9）過渡態(tài)如何驗(yàn)算虛頻呢？
比如一個6層原子層的slab上表面吸附小分子。slab下部3層原子是固定的。驗(yàn)算虛頻的時候，是不是還是固定下面三層原子，然后按照一般頻率計(jì)算方法來算虛頻？這樣的話，可以移動的原子數(shù)在20數(shù)量級上，考慮三個自由度，及其組合，就有很多很多可能了。請問該怎么設(shè)置這樣的過渡態(tài)虛頻計(jì)算呢？
16、關(guān)于概念的問題做個討論
（一）關(guān)于結(jié)合能。你說“結(jié)合能是定義為相距無窮遠(yuǎn)的原子結(jié)合形成一定結(jié)構(gòu)的物質(zhì)所放出的能量”
你和我說的沒區(qū)別，我說的是結(jié)合能是相對于“孤立原子做參考點(diǎn)的”，也就是它與周圍任何原子沒有相互作用，和你所說的相距無窮遠(yuǎn)一回事，我這個好像沒有任何錯誤。
===============
這里你說的是沒有錯誤，但是我覺得有必要先澄清一下。
（二）關(guān)于單點(diǎn)能。你說“它是第一性原理計(jì)算直接得到的能量，或者說是贗能，是一個空間點(diǎn)陣平均每陣點(diǎn)上采用贗勢計(jì)算所得到的能量，其中包含了結(jié)合能的貢獻(xiàn)，但是更多的，也包含了靠近芯區(qū)附近的電子在采用贗勢近似下的能量，這一部分能量既不是原子芯區(qū)附近電子能量的真實(shí)反應(yīng)，也不會影響化學(xué)鍵性質(zhì)，不會對結(jié)合能有所貢獻(xiàn)”。
我贊同你的大部分觀點(diǎn)，也提出你說的幾點(diǎn)錯誤，單點(diǎn)能準(zhǔn)確的來說它包含了所有哈密頓的量，而且這兒的單點(diǎn)能不是你所說的“平均每個原子的能量”，而是你計(jì)算的整個原胞的能量。但是這個能量有一個參考點(diǎn)。你可以看候博得，也可以看我回的下一個貼子，至于“影響不影響成鍵之類的內(nèi)容”固體力學(xué)上已經(jīng)說的很清楚了。
==============
從你的回復(fù)中，我可以知道你肯定沒有學(xué)過晶體學(xué)或者空間群理論，你應(yīng)該看看晶體學(xué)國際表中對于陣點(diǎn)的定義，陣點(diǎn)并不是每個原子，這里你的理解有問題，陣點(diǎn)是一個抽象點(diǎn)，一個晶體中包含所有對稱性的可以僅通過平移來構(gòu)造整個晶體的結(jié)構(gòu)所占據(jù)的位置就是一個陣點(diǎn)，換句話說，一個陣點(diǎn)，就是一個滿足平移對稱性的原子集團(tuán)，且該集團(tuán)內(nèi)部的位置滿足該晶體結(jié)構(gòu)的全部對稱性，而且它不僅僅是“原胞”
（三）你說“Free energy TOTEN是體系總能，要減去陣點(diǎn)上分布的原子的能量再除以平均原子數(shù)才是結(jié)合能（當(dāng)然，這個和你的計(jì)算腳本的設(shè)計(jì)有關(guān)），而且這還沒有考慮不加展寬時沒有被計(jì)算到的能帶的因素”
我不贊同你后面說的幾點(diǎn)。Free energy TOTEN從字面意思上我們也知道它的結(jié)果是自由能，你可以說它是總能，因?yàn)楦鶕?jù)我上面的推導(dǎo)，它們至少在數(shù)值是相等得。不在于你把它說成什么，你就是把它說成總能，其實(shí)它還是等于結(jié)合能，等于自由能，等于內(nèi)能。至于除不除原子總數(shù)在于你想得到的是平均每個原子的還是總體系的，這在于個人�？紤]不考慮展寬，那要看ISMEAR等于幾，做幾個實(shí)例就會感覺到它考慮沒考慮了。
===============
這個能量確切的說應(yīng)該是叫做考慮電子振動熵的體系總自由能，當(dāng)不考慮展寬的時候，它是等于總能的，如果你讀過Vasp的代碼，就知道TOTEN在vasp的計(jì)算中就是總能，這個和結(jié)合能不是一個概念，還包含有非成鍵部分的貢獻(xiàn)，至于內(nèi)能的定義，如果你閱讀過塞茲的現(xiàn)代固體理論，或者Pauling的書，或者讀過Morse當(dāng)時提出morse勢的那篇文獻(xiàn)，就應(yīng)該知道，固體物理中所使用的內(nèi)能，指的是離子實(shí)的動能和原子的“結(jié)合能”之和——這里結(jié)合能之所以要打引號，是因?yàn)榘炊x，是要形成穩(wěn)定結(jié)構(gòu)或者亞穩(wěn)態(tài)結(jié)構(gòu)時才能稱之為結(jié)合能，內(nèi)能定義并沒有考慮粒子芯區(qū)附近電子能量的影響，正如你所說，是“相對于“孤立原子做參考點(diǎn)的””，在0K下已經(jīng)不考慮動能，因此就應(yīng)是總能減去孤立原子的能量和才行，至于結(jié)合能，則是穩(wěn)定狀態(tài)下結(jié)構(gòu)的這個能量。
（四）你說“是否考慮展寬和結(jié)合能的定義沒有關(guān)系”
我也沒說和它的定義有什么關(guān)系啊，但是由于數(shù)學(xué)處理帶來的誤差，它對結(jié)合能的結(jié)果有一定的影響啊。
（五）對單個原子的結(jié)合能的計(jì)算應(yīng)該只計(jì)算Gamma點(diǎn)的能量，且用削除簡并。
你說得第五點(diǎn)我不懂，我算單個原子的能量時一直是按三個表面的構(gòu)造方法來算的，也沒想過什么簡并。還望有高手給我?guī)椭谖妩c(diǎn)怎么理解。
=============
簡單的操作是計(jì)算單個原子能量只考慮Gamma點(diǎn)，然后三邊都設(shè)置在10A以上，且不相等
至于原因，應(yīng)該去查量子力學(xué)的書，記得本科的時候，老師都會講到的。

根據(jù)F=U-TS，E（結(jié)合）=E（bulk）-nE(單個離散)。而E（bulk）就是U，通常我們?nèi)（離散）為參考點(diǎn)。也就是把E（離散）看為0，這樣推出來，在T=0時，F(xiàn)=E（結(jié)合）。它是包含你所說得那些，而且就像我在前面的貼子中說得，它就是總得H得到的能量，但是它有一個參考，而這個參考就是離散原子的能量
17、用vasp軟件研究表面小分子的吸附解離遇到幾個問題。
1 關(guān)于反應(yīng)物和產(chǎn)物的結(jié)構(gòu)，請問你們是如何構(gòu)建的？（參考文獻(xiàn)嗎？）能不能介紹一下相關(guān)的經(jīng)驗(yàn)，個人認(rèn)為好的開端是成功的一半，所以需要更多的經(jīng)驗(yàn)幫助。
答：.據(jù)說現(xiàn)在流行的過渡態(tài)算法多用CINEB，vasp自帶的NEB也可以，但是鏡像受力優(yōu)化和勢壘精確度以及收斂速度略遜于CINEB，lz可以google到其官方網(wǎng)站了解。初始和終了的結(jié)構(gòu)是局部能量最低的構(gòu)型（vasp做弛豫收斂的構(gòu)型），至于lz所關(guān)心的是從具體結(jié)構(gòu)之間的過渡的化，一是猜測，二是文獻(xiàn)吧�；境霭l(fā)點(diǎn)是能量趨向降低的構(gòu)型。
2 關(guān)于結(jié)構(gòu)的優(yōu)化，能量收斂和力收斂的標(biāo)準(zhǔn)是否應(yīng)該比默認(rèn)值高，有利于過渡態(tài)的搜尋？
答：DFT理論計(jì)算能量可能更準(zhǔn)確，對力的計(jì)算由于是對能量再求導(dǎo)，涉及數(shù)值算法原因可能不準(zhǔn)，所以收斂標(biāo)準(zhǔn)基本采用默認(rèn)即可，有時適當(dāng)提高（一個量級以內(nèi)）也是可以的。在合理的精度范圍內(nèi)討論出合理的結(jié)果就達(dá)到目的了，這是師兄告我的。基本上收斂精度高于默認(rèn)值去跑的話，那就god bless you了。
3 關(guān)于結(jié)構(gòu)虛頻的處理，通過學(xué)習(xí)了解了消除虛頻的方法--將對應(yīng)原子的坐標(biāo)加到原坐標(biāo)上。由于vasp的振動模式不能可視化，所以請問如何找到對應(yīng)的原子？另外關(guān)于虛頻，大家是如何處理的呢。
答：.虛頻沒有關(guān)注過。只知道在鞍點(diǎn)位置才應(yīng)該有一個虛頻，屬于鞍點(diǎn)具有的特性吧。為何要消除呢？根據(jù)這個好像才能按速率理論算出反應(yīng)率的吧，這個我沒有經(jīng)驗(yàn)。只關(guān)注過勢壘。
4 關(guān)于images的問題，最近利用腳本產(chǎn)生了2、4的POSCAR的文件，發(fā)現(xiàn)只有2的能跑起來，我的服務(wù)器是4個核的�？磩e人發(fā)的貼：指定4個cpu同時計(jì)算，應(yīng)該是1個cpu一個image 。為什么我設(shè)置4個image跑不起來呢？
答：image數(shù)量越多，相互之間受力需要同時滿足收斂準(zhǔn)則，速度肯定會降低，但是這樣找到的能量最低路徑也會與實(shí)際更加符合。沒錯，每個image事實(shí)上是需要相當(dāng)于單獨(dú)計(jì)算的構(gòu)型的cpu的數(shù)目的，所以4個核的話，平均一個image一個cpu，實(shí)在是不一定跑得動。如果體系原子數(shù)較多，k點(diǎn)網(wǎng)格較密，能量截?cái)嘀翟僖桓撸烙?jì)經(jīng)過很久才可能有結(jié)果。
5、在計(jì)算dos時，為什么會出現(xiàn)WARNING:stress and forces are not correct。這是那些原因造成的，計(jì)算的是多鐵物質(zhì)。還有就是在優(yōu)化時，TOTAL-FORCE這項(xiàng)怎么全是0啊，原子這個沒有關(guān)系的
答：用優(yōu)化后的CONTCAR算band或dos都會出現(xiàn)這個提示位置沒有變化啊，這是怎么回事��？
6、我在做能帶計(jì)算，發(fā)現(xiàn)自洽計(jì)算和非自洽計(jì)算中，OUTCAR文件中都可以找到一個費(fèi)米能級的大小，但是現(xiàn)在迷茫的是，不知道在畫能帶時減去的費(fèi)米能級的大小應(yīng)該取哪個結(jié)果文件中得值？是取自洽計(jì)算的結(jié)果文件中的E-fermi ，還是取非自洽的結(jié)果文件中的E-fermi 的大��？？
答：算能帶的那次計(jì)算，指的是自洽（求解WAVECAR和CHGCAR的那一步）的那步計(jì)算，還是非自洽（直接讀取自洽的電荷和波函數(shù)，進(jìn)行非自洽計(jì)算）的那次計(jì)算啊？我對這個說法的標(biāo)準(zhǔn)不是很理解，希望您能再指導(dǎo)一下，我覺得畫能帶圖時，費(fèi)米能級的取值是否正確是很關(guān)鍵的，取自洽計(jì)算得到的E-fermi

3.第一性原理計(jì)算Boltztrap篇，里面也有大量VASP介紹
（http://www5.hp-ez.com/hp/calculations/page104）
4.vasp學(xué)習(xí)資料參考
http://www.gaoyang168.com/bbs/viewthread.php?tid=7474353&fpage=1
http://www.gaoyang168.com/bbs/viewthread.php?tid=596999&fpage=1
以上內(nèi)容僅供參考
不要說話
50樓2014-11-21 08:54:33
已閱   回復(fù)此樓   關(guān)注TA 給TA發(fā)消息 送TA紅花 TA的回帖
相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 ljw4010 的主題更新
普通表情 龍 兔 虎 貓 高級回復(fù) (可上傳附件)
最具人氣熱帖推薦 [查看全部] 作者 回/看 最后發(fā)表
[考研] 336求調(diào)劑 +4 收到VS 2026-03-20 4/200 2026-03-23 19:02 by macy2011
[考研] 上海電力大學(xué)材料防護(hù)與新材料重點(diǎn)實(shí)驗(yàn)室招收調(diào)劑研究生（材料、化學(xué)、電化學(xué)，環(huán)境） +3 我愛學(xué)電池 2026-03-23 3/150 2026-03-23 17:16 by AZMK
[考研] 一志愿南京理工大學(xué)085701資源與環(huán)境302分求調(diào)劑 +5 葵梓衛(wèi)隊(duì) 2026-03-18 7/350 2026-03-23 16:26 by lingjue
[考研] 石河子大學(xué)（211、雙一流）碩博研究生長期招生公告 +3 李子目 2026-03-22 3/150 2026-03-22 21:01 by 怎么釋懷
[考研] 319求調(diào)劑 +4 小力氣珂珂 2026-03-20 4/200 2026-03-22 15:53 by ColorlessPI
[考研] 求調(diào)劑 +5 Zhangbod 2026-03-21 7/350 2026-03-22 13:13 by Zhangbod
[考研] 286求調(diào)劑 +10 Faune 2026-03-21 10/500 2026-03-21 23:34 by 314126402
[考研] 求調(diào)劑 +3 13341 2026-03-20 3/150 2026-03-21 18:28 by 學(xué)員8dgXkO
[考研] 0703化學(xué)297求調(diào)劑 +3 Daisy☆ 2026-03-20 3/150 2026-03-21 17:45 by ColorlessPI
[考研] 268求調(diào)劑 +9 簡單點(diǎn)0 2026-03-17 9/450 2026-03-21 15:37 by lature00
[考研] 070300化學(xué)319求調(diào)劑 +7 錦鯉0909 2026-03-17 7/350 2026-03-21 03:46 by JourneyLucky
[考研] 二本跨考鄭大材料306英一數(shù)二 +3 z1z2z3879 2026-03-17 3/150 2026-03-21 02:29 by JourneyLucky
[考研] 324分 085600材料化工求調(diào)劑 +4 llllkkkhh 2026-03-18 4/200 2026-03-21 01:24 by JourneyLucky
[考研] 一志愿華中科技大學(xué)，080502，354分求調(diào)劑 +5 守候夕陽CF 2026-03-18 5/250 2026-03-21 01:06 by JourneyLucky
[考研] 一志愿 西北大學(xué) ，070300化學(xué)學(xué)碩，總分287，雙非一本，求調(diào)劑。 +3 晨昏線與星海 2026-03-18 3/150 2026-03-21 00:46 by JourneyLucky
[考研] 296求調(diào)劑 +6 www_q 2026-03-18 10/500 2026-03-20 23:56 by JourneyLucky
[考研] 一志愿中海洋材料工程專碩330分求調(diào)劑 +8 小材化本科 2026-03-18 8/400 2026-03-20 23:16 by JourneyLucky
[考研] 330求調(diào)劑 +4 小材化本科 2026-03-18 4/200 2026-03-20 23:13 by JourneyLucky
[考研] 招收調(diào)劑碩士 +4 lidianxing 2026-03-19 12/600 2026-03-20 12:25 by lidianxing
[考研] 320求調(diào)劑0856 +3 不想起名字112 2026-03-19 3/150 2026-03-19 22:53 by 學(xué)員8dgXkO
信息提示
請?zhí)钐幚硪庖?/div>
欧美大鸡吧男操女啊啊啊视频 | 99久久免费播放在线观看视频| 天堂网免费在线电影| 欧美国产精品久久久免费| 天海翼亚洲一区在线观看| 正在播放麻豆精品一区二区| 伊人久久综合国产精品| 亚洲gay视频在线观看| 一级做性色a爱片久久片| 亚洲欧洲无码一区2区无码| 青娱乐免费最新视频| 日本清纯中文字幕版| 最新激情中文字幕视频| 69精品人妻久久久久久久久久久| 伊人久久综合国产精品| 一区二区三区高清视频3| 2020国产激情视频在线观看| 久久久久高潮白浆久久| 五月的婷婷综合视频| 男人和女人的逼视频| 国产精品乱码高清在线观看h| 国产精品网站的黄色| 538欧美在线观看一区二区三区 | 亚洲综合在线视频在线播放| 午夜福利午夜福利影院| 日日躁夜夜躁狠狠操| 久久人人爽人人爽人人av东京热| 国产资源网站在线播放| 人妻女侠被擒受辱记| 精品不卡一区二区三区| 天天干天天操天天日天天日| 又粗又长又硬又黄又爽| 国产免费久久精品99re丫丫| 91精品麻豆91夜夜骚| 亚洲国内精品久久久久久久| 国产原创一区二区三区在线播放| 一二三四区国产在线观看| 精品人妻在线激情视频| 裸露视频免费在线观看| 日本东京热视频欧美视频| 亚洲第一成年偷拍视频| 鸡巴在里面福利视频在线观看| 久久亚洲国产成人精品麻豆| av无限看熟女人妻另类av| 欧美日本在线免费视频| 日本高清在线观看不卡视频| 熟妇高潮久久久久久久| 99久久人人爽亚洲精品美女| 亚洲全国精品女人久久久| 91精品国产人妻麻豆| 在线免费观看a视频免费| 港台美女明星av天堂| 猫咪亚洲中文在线中文字幕| 日本老女人日比视频| 天天操天天干天天谢| 第一福利视频在线观看| —区二区三区女厕偷拍| 中文字幕av特黄毛片| 夜色福利视频免费观看| 大奶熟妇激情操逼逼| 国产经典精品欧美日韩| 国产一区二区三区四区精| 国语对白性爱三级片免费看| 亚洲av综合av一去二区三区| 中文字幕久久久国产| 亚洲制服丝袜美腿在线| 中国特黄色性生活片| 亚洲美女午夜激情视频在线观看| 18福利视频在线观看| 四季av人妻一区二区三区| 国产成人在线观看hd| 亚洲一区二区在线视频观看免费| 中字幕人妻熟女人妻a62v网| 亚洲午夜高清在线观看| 成人精品动漫一区二区| 夜夜人人干人人爱人人操| 亚洲欧洲一区二区三区在线| 午夜精品一区二区三区不卡顿 | 无人区一码二码三码区别在哪| 成人十欧美亚洲综合在线 | 亚洲av中文无码网站| 啊不行啊操逼好爽大鸡吧视频| 中文字幕在线免费观看成人| 自拍丝袜国产欧美日韩| 成人黄色录像在线观看| 亚洲男人天堂最新网址大全| 青青在线视频看看| 亚洲欧美一级特黄大片| 97香蕉久久国产超碰| 青青操天堂在线观看视频| 韩国毛片w妈妈的朋友7| 正在播放麻豆精品一区二区| 国产精品剧情av在线播放| 蜜乳av一区二区三区免费观看| 夜夜骚av一二三区| 亚洲av毛片一区二区三区网| 男女69视频在线观看免费| 天天在线播放日韩av| 久久99国产中文丝袜| 天天碰天天摸天天搞| 日本东京热视频欧美视频| 91超碰九色porny| 久久中文字幕av一区二区| 黄版视频在线免费观看| 伊人情人成综合视频| 69视频在线精品国自产拍| 女人扒开逼让男人操| 成人超碰一区二区三区| 中文字幕观看中文字幕免费 | 啪啪啪网站免费在线看 | 青娱乐免费最新视频| 欧美亚洲愉拍一区二区三区| 国产主播诱惑毛片av| 一区二区三区高清视频3| 92午夜免费福利视频www| 欧美日韩亚洲国产视频二区| 大屁股熟女一区二区视频| 啪啪啪网站免费看视频| 男女爱爱好爽视频免费看| 亚洲制服丝袜美腿在线| 99精品久久精品一区二区| 国产伦理二区三区在干嘛呢| 全彩漫画口工18禁| 99999久久久精品| 在宿舍强奷两个清纯校花| 欧美日本在线免费视频| 女生抠逼自慰啊啊啊啊啊啊啊下载| 亚洲成人五月婷婷久久综合| 97人妻av人人澡人人爽| 中文字幕久久久国产| 亚洲乱码国产乱码精品精视频| 亚洲成人av在线一区二区| 99福利一区二区视频| 日韩成人精品久久久免费看| 69国产在线视频网站| 亚洲免费在线不卡视频| 熟女国内精品一区二区三区| 自拍丝袜国产欧美日韩| 女同大尺度视频网站在线观看| 91 精品视频在线看| 天天夜夜久久精品综合| 宅男噜噜噜66国产在线观看| 国产高清自拍偷拍在线| 亚洲乱码国产乱码精品精视频| 成人十欧美亚洲综合在线| 美女欧美视频在线观看免费| 亚洲成人激情在线综合| 黑人大巨屌操美女逼| 玖辛奈18禁同人污本子| 美国伦理片午夜理论片| 东京热男人的天堂视频| 久久久精品人妻无码专区不卡| 熟妇人妻丰满久久久久久久| av中文字幕国产精品| 69av精品国产探花| 在线观看2022av| 制服丝袜中文字幕熟女人妻| 懂色av之国产精品| 不卡高清一区二区三区| 亚洲|久久久久久一二三区丝袜| 天天插天天操天天射天天干| 日本午夜福利免费在线播放| 天天操天天搞天天操| 日韩av水蜜桃一区二区三区| 亚洲成人欧洲成人在线| 在线成人教育平台排名| 成人超碰一区二区三区| 亚洲国产精品青青草| 欧美成人红桃视频在线观看| 2026天天操天天干| 亚洲中文字幕在线视频观看二区| 91九色pony蝌蚪| 中文字幕 首页 人妻| 97精品视频,全部免费| 欧美亚洲国产一区二区| 亚洲激情噜噜噜久久久| 午夜宅男电影av网站| 亚洲AV无码久久精品国产一区老| 老司机免费视频福利0| 亚洲蜜桃久久久久久| 2026天天操天天干| 亚洲熟妇丰满多毛xxxx网站| 欧美一区二区三区爽爽| —区二区三区女厕偷拍| 欧美成人短视频在线播放| 国内精品一区二区2021在线| 欧美精品999不卡| 亚洲精品国产99999| 女女抠逼白虎白丝袜| 亚洲全国精品女人久久久| 熟妇人妻丰满久久久久久久| 免费在线小视频你懂的| 欧美色视频网址大全| 欧美一级日韩一级亚洲一级va | 亚洲午夜高清在线观看| 超级黄肉动漫在线观看| 国产精品 亚洲欧美 自拍偷拍| 一区二区三区免费版在线| 成人午夜麻豆大胆视频| 999精品视频免费在线观看| 五十岁熟妇高潮喷水| 国产资源在线观看二区| 亚洲欧美精品日韩偷拍| 日韩av熟妇在线观看| 中文字幕日本一二三区| 日本成年视频在线免费观看| 日韩黄色在线观看网站上| 久草久热这里只有精品| 青青在线视频看看| 美女激情久久久久久久| 亚洲国产综合久久精品| 奇米网首页神马久久| 呻吟求饶的人妻中文字幕| 亚洲中文字幕最新地址| 国产av在线免费视频| 天堂av在线最新地址| 5566熟女人妻人妻| 婷婷色综合五月天视频| 天天日天天干天天日天天干天天| 中文字幕欧美一区二区视频| 婷婷色综合五月天视频| 国产农村乱子伦精精品视频| 抽插小穴啊啊啊视频| 无码精品黑人一区二区老人| 黄在线看片免费人成视频| 亚洲|久久久久久一二三区丝袜| 欧美在线观看一区二区不卡| 区一区二区三免费观看视频| 不卡一区二区视频在线| 国产91精品福利系列| 在线观看2022av| 91系列视频在线播放| av大尺度一区二区三区| 大奶熟妇激情操逼逼| 91精品麻豆91夜夜骚| 自拍偷自拍亚洲精品10p| 38av一区二区三区| 久久久久久久久久久久久国产| 精品欧美乱码久久久| 国产三级自拍视频在线观看网站| 亚洲精品激情视频在线观看| 成人做爰av在线观看网站| 亚洲综合第一区二区| 福利美女视频在线观看| 91porny九色视频偷拍| 核xp工厂精品久久亚洲| 99 re国产精品| 日本高清激情乱一区二区三区| 日韩成人精品久久久免费看| 全国熟妇精品一区二区免费视频| 桃色成人开心激情网| 亚洲天堂男人的天堂| 可在线免费观看av| 人妻激情偷乱一区二区三区av| 国产漂亮白嫩美女在线图片 | 国产资源网站在线播放| 日本人妻熟妇丰满成熟HD系列 | 黑人大吊大战亚洲女人。| 午夜久久久久欠久久久久| 高清av在线婷一区二区色日韩| 首页欧美日韩中文字幕| 国产av高清二区三区| 青娱乐免费视频一二三| 日本男女免费福利视频| 成人大片男人的天堂| 亚洲一区二区偷拍女厕所| 日本一区二区高清av中文| 极品少妇高潮喷水日出白浆| 久久久人妻免费视频| 户外露出视频在线观看| 91色哟哟视频在线观看| 制服丝袜 中文字幕 日韩| 77亚洲视频在线观看| 亚洲天堂男人的天堂| 欧美久久蜜臀蜜桃资源吧| 日本四十路人妻熟女| 岳的大肥屁熟妇五十路| 亚洲avav天堂av在线网毛片| 久久99精品热在线观看| 538欧美在线观看一区二区三区 | 亚洲三级综合在线观看| 国产中年夫妇激情高潮| 一区二区三区av免费天天看| 国产欧美福利在线观看| 天天看片天天摸天天操| 日本午夜福利免费在线播放| 日韩av水蜜桃一区二区三区| 最新免费在线观看污视频| 日本清纯中文字幕版| 91精品国产91久久久久久密臀| 77亚洲视频在线观看| 天天操天天干天天谢| 男女插鸡巴视频软件| 欧美日韩久久丝袜在线| 成人午夜高清福利视频| 大乳丰满人妻中文字幕韩国hd| 日韩激情文学在线视频| 亚洲欧美成人激情在线| 中文字幕人妻一区色偷偷久久| 182tv精品免费在线观看| 国产激情免费在线视频| 交换的一天中文字幕在线视频| jizzjizz国产精品传媒| 丰满放荡熟妇在线播放| 极品风骚人妻3p视频| 亚洲av手机免费在线| lutu玩弄人妻短视频| 国产视频成人自拍蝌蚪视频| 亚洲av综合av一去二区三区| 97精品国产91久久久| 午夜92福利1000| 久久一级片三上悠亚| 午夜夫妻性生活视频| 亚洲欧美精品日韩偷拍| 亚洲自拍偷拍一区二区中文字幕| 色网站在线观看免费| 欧美日本国产一区二区| 黑吊操欧美极品美女| 男人和女人的逼视频| 午夜精品久久久久久久精品乱码| 日产国产欧美精品另类| 久久99久久99久久97的人| 91麻豆精品国产在线| av毛片在线观看网址| 日韩一区二区在线播放观看| 日本欧美视频在线免费| 日本欧美国产在线一区| 美女黄色啊啊啊啊视频| 中文字幕人妻一区色偷偷久久 | 中文字幕日韩人妻在线三区 | 久草视频在线看免费| aaaa级少妇高潮在线观看| 顶级欧美色妇xxxx| jiee日本美女视频网站| 夜夜操天天干夜夜操| 国产男人的天堂一区| 亚洲一区亚洲二区成人福利| 中文字幕人妻一区二区视频系列 | 999国产精品视频免费看| 制服丝袜 中文字幕 日韩| 日韩av熟妇在线观看| 中文字幕在线免费观看人妻 | 亚洲免费在线不卡视频| 中文字幕熟女人妻丝袜丝在线| 在线能看视频你懂的| 亚洲成年人精品国产| 亚洲 综合 欧美 一区| 中日韩又粗又硬又大精品| 看女人大BB群伦交| 最新国产精品综合网高清| 亚洲avav天堂av在线网毛片| 中文字幕熟女人妻一区| 在线观看中文字幕视频成人| 日本东京热最新中文字幕| 户外露出视频在线观看| 韩国一级片最火爆中文字幕| 38av一区二区三区| 美女一区二区四区六区八区| 大秀成年人国产精品视频| 69久久夜色精品国产69乱电影| 黑人大巨屌操美女逼| 一二区二区不卡视频| 成人av在线视频免费| 玖玖资源站在线观看亚洲| 日本一道中文字幕99| 女人的天堂 av在线| 日本成年视频在线免费观看| 欧洲亚洲一区二区三区四区| 18福利视频在线观看| 中文字幕福利视频在线一区| 69精品互换人妻4p| 欧美日韩不卡视频合集| 天天色天天射天天日天天干| 91大神福利视频网| 亚洲国产精品久久久久久无码| 成年人免费福利在线| 高清av在线婷一区二区色日韩| 亚洲最大先锋资源采集站| 黑人黄色免费一级av| 全国熟妇精品一区二区免费视频 | 内地精品毛片在线观看| 插鸡视频免费网站在线播放| 夜夜躁av麻豆男| 68视频在线免费观看| jiee日本美女视频网站| 午夜偷拍的视频久久久免费大全| 午夜宅男电影av网站| 中文字幕福利视频第四页| 国产不卡免费在线观看| 欧美一区二区三区爽爽| 欧美日本在线免费视频| 亚洲色图日韩在线视频观看| 亚洲国产精品一区二区第二页| 天天日天天亲天天操| 999国产精品视频免费看| www国产亚洲精品久久久| 一区二区三区高清视频3| 新亚洲天堂男子av| 一区二区三区国产精华液区别大吗| 国产男女无套?免费网站下载 | 黑人爆操女人免费视频| 日本成年视频在线免费观看| 人妻少妇的va视频| 亚洲精品国品乱码久久久久| av在线中文字幕在线| 日韩成人精品久久久免费看| 国产探花自拍亚洲av| 亚洲综合首页综合在线观看 | 91精品国产人妻麻豆| 国产精品成人免费电影| 亚洲成人,国产精品| 国产igao激情在线视频入口| 国内销魂老女人老泬| 高清av在线婷一区二区色日韩| 欧美日韩亚洲国产视频二区| 97人妻在线视频自拍| 亚洲人成大片在线观看| 插鸡视频免费网站在线播放| 蜜桃臀少妇白色紧身裤细高跟| 日韩在线 中文字幕| 日韩av熟妇在线观看| 中文字幕 首页 人妻| 亚洲国产精品青青草| 免费看一级高潮喷水片| 核xp工厂精品久久亚洲 | 中文字幕观看中文字幕免费 | 97cao在线视频| 国产最新av在线免费观看| 亚洲精品乱码久久久久app| 最新日韩av电影在线播放| 青娱乐免费视频一二三| 亚洲最强的25个城市| 日本特级黄片免费观看| 欧美一级日韩一级亚洲一级va| 99re这里是国产精品首页| 精品av天堂毛片久久久| 中文人妻av一区二区三区| 97视频人人爱麻豆| 成人人妻h在线观看| av天堂hezyo| 超级黄肉动漫在线观看| 不卡高清一区二区三区| 国产毛片特级Av片| 91精品麻豆91夜夜骚| 亚洲一区二区在线激情| 亚洲成年人精品国产| 久久精品国产亚洲av清纯| 国产剧情av在线免费观看| 午夜美女福利视频在线| 精品av天堂毛片久久久| 色欲AV蜜桃一区二区三| 美女把腿张开给男的捅| 久久久久夜色国产精品电影| 婷婷色九月综合激情丁香| 日韩国产欧美一区二区三区粉嫩| 污网址在线观看视频| 4438全国成人免费视频| 久久综合狠狠综合久久综| 欧美肥妇久久久久久| 91人妻人人爽色啊啊啊| 成人人妻h在线观看| 男人av一区二区三区| 欧美成人屋影院在线视频观看| 91福利高清在线播放| 91超碰九色porny| 亚洲中文字幕在线视频观看二区| 一级做性色a爱片久久片| 偷拍欧美日韩另类图片| 有码一区二区三区四区五区| 亚洲妹妹我爱你在线观看 | 大乳丰满人妻中文字幕韩国hd| 欧美情色av在线观看| 天天摸天天干夜夜操| 77亚洲视频在线观看| 18在线观看免费观看| 69视频在线精品国自产拍| 国产福利一区二区三区在线观看 | 色狠狠色综合久久久绯色| 午夜精品秘一区二区三区| 性感人妻 中文字幕| 美女av色播在线播放| 成人精品影视一区二区| 天天操天天干天天谢| 亚洲少妇色小说综合| 日本一区二区三区区别| 在线成人教育平台排名| 日韩av电影中文在线免费观看| 国产伦理二区三区在干嘛呢| 日本老熟老熟妇七十路| 中文字幕欧美人妻在线.| 美女福利网站在线播放| 成人精品影视一区二区| 老司机免费视频福利0| 美女露阴道让男人捅| 天天摸天天干夜夜操| 国产精美视频精品视频精品| 夏目彩春av在线看| 久久久久久a女人处女| 51vv精品视频在线观看| 人人妻人人爽人人摸| 国产白丝一区二区三区av| 无码精品黑人一区二区老人| 99久久99九九九99九| 亚洲制服丝袜资源网| 亚洲av综合av一去二区三区| 欧美日韩亚洲tv不卡久久| 亚洲精品激情视频在线观看| 免费中文三级在线观看| av网页免费在线观看| 国产视频成人自拍蝌蚪视频| 午夜久久人妻一级内射av网址| 国长拍拍视频免费孕妇| 国产青青青青草免费在线视频| 啪啪啪网站免费看视频| 欧美男女一区二区三区| 中国特黄色性生活片| 婷婷六月天在线视频| 亚洲欧美综合另类最新| 抽插小穴啊啊啊视频| 欧美第一激情综合网欧美激情| 中文字字幕在线精品乱码| 亚洲无码专区中文字幕专区| 岛国av成人午夜高清| 国产91九色视频在线观看| 在线能看视频你懂的| 国产av精品一区二区三区久久| 国产精品蝌蚪自拍视频| 91香蕉国产亚洲一二三区| 夜夜操天天干夜夜操| 欧美一级aaaaaaa片| 欧美一区二区三区爽爽| 一区二区三区四区视频精品免费| 干逼又爽又黄又免费的视频| 欧美黑人1区2区3区| 夜夜爽夜夜操夜夜爱| 婷婷六月天在线视频| 亚洲宅男噜噜噜66在线观看| 国产91免费在线观看| 久草视频在线视频在线视频| 欧美操大黑鸡巴视频在线观看| 欧美亚洲愉拍一区二区三区| 一看就是假奶的av| 天天操天天干天天谢| 最新日韩av电影在线播放| 色哟哟亚洲乱码国产乱码精品精 | 青青草成人免费自拍视频| 人妻系列中文字幕大乳丰满人妻 | 亚洲avav天堂av在线网毛片| 亚洲欧美精品海量播放| 亚洲第一区av中文字幕| 50熟妇一区二区三区| 国产做A爱免费视频在线观看| 自拍偷拍亚洲综合第一页| 一区二区三区婷婷中文字幕| 青青操久久综合激情| 国产高清视频www夜色资源| 欧美日本亚欧在线观看| 超碰在线免费观看视频97| 欧美亚洲国产一区二区| 丰满放荡熟妇在线播放| 性感人妻 中文字幕| 一区二区三区 国产日韩欧美| 久久中文字幕av一区二区| 婷婷综合缴情亚洲五月伊人| 国产伦理二区三区在干嘛呢| 亚洲avav天堂av在线网毛片| 欧美视频亚洲视频在线| 天天干天天日天天弄| 国产男人的天堂一区| 成人资源中文在线观看| 亚洲综合第一区二区| 亚洲|久久久久久一二三区丝袜| 久久人妻诱惑我视频| 午夜久久久久欠久久久久| 大尺度av毛片在线网址| 久久久视频在线播放| 综合久久伊人久久88| 99re这里是国产精品首页| 偷拍熟女大胆免费视频| 鸡巴在里面福利视频在线观看| yellow在线亚洲精品一区| 午夜国产免费视频亚洲| 男生用大肌巴操美女骚穴| 国产精品网站亚洲发布| 真人一进一出抽搐大尺度视频| 熟妇人妻av无码中文字幕| 国产资源在线观看二区| 可以免费观看日韩av| 亚洲成人激情在线综合| 欧美操大黑鸡巴视频在线观看| 精品国产污污污免费入口| 人妻色综合aaaaaa网| 国产午夜羞羞一区二区三区| 视频自拍偷拍视频自拍| 天天干夜夜操91视频网站| julia人妻av一区二区三区| 欧美日韩一区二区三区成人影院| 天天干夜夜操91视频网站| 日本高清激情乱一区二区三区| 亚州av嫩草av极品在线观看| 午夜精品久久久久久久精品乱码| 久久久久国产精品二区| 91麻豆精品国产在线| 欧美男男在线观看视频网站| 中文字幕观看中文字幕免费| 成人av中文字幕在线看| 亚洲精品9999蜜桃| 欧美巨大另类极品video| 国产白丝一区二区三区av| www国产亚洲精品久久久| 国产av精品一区二区三区久久| 国产精品内射婷婷一级| 国产精品成人免费电影| v天堂国产精品久久| 天天操天天舔天天爽| 青青草成人免费自拍视频| 亚洲精品1卡2卡3卡| 少妇精品视频一区二区免费看| 国产 亚洲 欧美 自拍| 天天干天天色综合久久| 欧美成人少妇人妻精品| 天天搞天天操天天干| 午夜情色一区二区三区| 亚洲成人av在线一区二区| 蜜臀久久精品久久久久久av | 亚洲av三级电影在线观看| 91国产精品乱码久久久久久| 欧美精品一区二区三区观看| 99久久国产精品免费热| 午夜美女福利视频在线| 免费看日韩黄视频在线观看| 中文字幕av特黄毛片| 免费啪啪啪网站在线观看| 农村大炕有肉大屁股熟妇| 男生和女生羞羞91在线看| 欧美一级日韩一级亚洲一级va| 成人av在线视频免费| 五十岁熟女高潮喷水| 日韩最近中文在线观看| 日本欧美视频在线免费| 在宿舍强奷两个清纯校花| 天天干夜夜操91视频网站| 亚洲熟女一区二区六区| 精品视频在线观看免费99| 上床啪啪啪免费视频| 亚洲无码专区中文字幕专区| 大屁股熟女一区二区视频| 狠狠操狠狠操狠狠插| av福利免费体验观看| 天天操天天日天天碰| 日韩人妻中文字幕二区| jiee日本美女视频网站| 二十四小时日本高清在线观看| 久久一级片三上悠亚| 欧美成人少妇人妻精品| 欧美精品乱码99久久蜜桃免费| 伊人久久综合国产精品| 欧美亚洲国产一区二区| 国产av啊啊啊啊啊啊啊| 成人大片男人的天堂| 国色天香一二三期区别大象| 黑鸡巴肏少妇逼视频| aaaa级少妇高潮在线观看 | 夜色福利视频免费观看| 成人人妻h在线观看| 青青在线视频看看| 亚洲av 综合av| 伦理在线观看未删减中文字幕| 熟妇精品午夜久久久久| 黄色大片一级老太太操逼| 日韩人妻中文字幕区| 快使劲弄我视频在线播放| 最近最新最好看的中文字幕| 午夜福利国产精品久久久久| 精品久久久久久久久久久久久| 亚洲自拍偷拍一区二区中文字幕| 亚洲熟女人妻自拍在线视频| 裸日本资源在线午夜| 国产成人在线观看视频播放| 99国产精品久久99久久久| 欧美区一区二区三视频| 亚洲一区二区三区国产精品电影| 婷婷一区二区三区五月丁| 老鸭窝在线毛片观看免费播放 | 欧美日韩福利视频网| 青娱乐这里只有精品| 国产乱码有码一区二区三区| 亚洲一区二区三区国产精品电影 | 一区二区三区国产在线成人av| 日本国产亚洲欧美色综合| 亚洲成人五月婷婷久久综合| 性感美女人妻久久久| 男女啪啪啪啪91av日韩| 亚洲成人五月婷婷久久综合| 亚洲一区视频中文字幕在线播放| 久久人人爽人人爽人人av东京热| —区二区三区女厕偷拍| 亚洲国产精品自拍偷拍视频在线 | 91九色pony蝌蚪| 亚洲国内精品久久久久久久| 综合激情网,激情五月| 中文字幕福利视频在线一区| 亚洲国产精品自拍偷拍视频在线| 免费的啪啪视频软件| 手机视频在线观看一区| 69精品互换人妻4p| 涩涩黄片在线免费观看| 视频免费在线观看网站| 午夜国产一区二区三区| 久久一级片三上悠亚| 91香蕉国产亚洲一二三区| 4日日夜夜精品视频免费| 天天操天天舔天天射天天日天天干 | 女同性恋av在线播放| 68福利精品在线视频| 人妻中文字幕亚洲在线| 亚洲精品乱码久久久久app| 亚洲一区二区精品三区视频| 亚av一二三在线观看| 自拍偷拍 国产激情| 中国特黄色性生活片| 欧美国产精品久久久免费| 午夜精品视频免费观看| 999国产精品视频免费看| 亚洲午夜熟女在线观看| 成人黄色录像在线观看| 日本男女免费福利视频| 丝袜美腿日韩av一区| 亚洲欧美成人激情在线| 国际精品熟女一区二区| 熟女阿高潮合集一区二区| 国产主播诱惑毛片av| 日本有码精品一区二区三区| 日韩免费黄色片在线观看| 亚洲色图日韩在线视频观看| 欧美性感美女热舞视频| 韩国在线播放一区二区三区| 国产精品中文字幕丝袜| 亚洲另类欧美综合久久| 亚洲午夜精品视频节目| 东京热日韩av在线| 伊人情人成综合视频| 老司机免费视频福利0| 美女网站视频久久精品| 加勒比东京热绿帽人妻多人操| 亚洲欧洲一区二区三区在线| 亚洲av综合av一去二区三区| 快使劲弄我视频在线播放 | 污网址在线观看视频| 裸露视频免费在线观看| 97精品国产91久久久| 成人十欧美亚洲综合在线| 美女张开腿给男人桶爽的软件| 嗯~嗯~啊啊啊~高潮了软件| 男女插鸡巴视频软件| yellow在线亚洲精品一区| 亚洲国产电影的一区| 荣立三等功退休有什么待遇| 午夜8050免费小说| 国产原创一区二区三区在线播放| 核xp工厂精品久久亚洲| 无码精品黑人一区二区老人| 欧美aaaa性bbbbaaaa| 强乱人妻中文字幕日本| 免费绝清毛片a在线播放| 99re这里是国产精品首页| 92在线播放观看视频| 美女把腿张开给男的捅| 久久一级片三上悠亚| 最新中文字幕久久久久| 日本欧美国产在线一区| 快进来插我的逼嗯啊视频| 亚洲精品1卡2卡3卡| 核xp工厂精品久久亚洲| 91九色国产在线视频| 18福利视频在线观看| 欧美日本在线免费视频| 国内销魂老女人老泬| 少妇被粗大的猛进69视频| 久久中文字幕av一区二区 | 在线看日韩av不卡| 在线视频国产精品欧美| 人妻系列在线免费视频| 国产美女视频带a∨黄色片| 伊人精品成人综合网| 二十四小时日本高清在线观看 | 自拍偷自拍亚洲精品10p| 一区二区三区高清视频3| 杜达雄啪啪毛片视频| 亚洲精品国品乱码久久久久| 天天色 天天操 天天好逼| 国产天堂av不卡网| 天堂av国产av伦理av| 亚洲 偷拍 自拍 欧美| 人人妻人人澡人人爽97| 蜜乳视频一区二区三区| 亚洲一区二区在线视频观看免费 | 亚洲欧洲无码一区2区无码| 天天操天天干加勒比久久| 日本少妇精品免费视频| 91超碰九色porny| 老牛影视在线一区二区三区| xxnxx国产美女| 亚洲理论在线a中文字幕97| 中文字幕日韩人妻在线三区| 青青免费观看视频| 新亚洲天堂男子av| 精品久久久久久久久久久久久| 中文字幕熟女人妻一区| 九九视频在线观看全部| 欧美日韩亚洲国产视频二区| 日本不卡 中文字幕| 国产经典精品欧美日韩| 日本电影一级人妻在线播放四区| 松本菜奈实最新av在线| 91精品国产欧美在线| 国产做A爱免费视频在线观看| 一区二区三区四区影片| 国产精品网站亚洲发布| 欧美精品999不卡| 999久久久人妻精品一区| 欧美视频免费观看777| 亚洲国产日韩欧美一区二区三区,| av在线播放观看h| 中文字幕熟女人妻丝袜丝在线| 国产探花自拍亚洲av| 欧美日韩久久丝袜在线| 国产福利三级在线观看| 乌克兰美女操逼高清内射视频| 在线观看网站伊人网| 91麻豆精品国产在线| 欧美成人久久久桃色aa| 亚洲精品激情视频在线观看 | 国产精品乱码高清在线观看h| 欧美亚洲愉拍一区二区三区| 裸日本资源在线午夜| 美女激情久久久久久久| 首页欧美日韩中文字幕| 1区3区4区产品乱入视频| 午夜精品一区二区三区不卡顿| 一二三四区国产在线观看 | 91在线九色porny| 精品视频一区二区三区◇| 中文字幕熟女人妻丝袜丝在线| 日韩国产欧美久久一区| 91 精品视频在线看| 强乱人妻中文字幕日本| 欧美日韩综合精品无人区| 男人电影天堂在线观看| 欧美日本国产一区二区| 精品欧美乱码久久久| 中文字幕熟女乱一区二区| 亚洲av 综合av| 亚洲全国精品女人久久久| 最新中文字幕久久久久| 国产人妻熟女ⅹxx丝袜| 成人午夜av电影网| 欧美黄色性视频网站| 猫咪亚洲中文在线中文字幕| 后入日韩翘臀蜜桃臀美女| 99久久碰碰人妻国产| av里面的动作是真进去吗| 久久午夜免费鲁丝片| 一看就是假奶的av| 91色老久久精品偷偷蜜臀| 国模伊人久久精品一区二区三区| av男人站在线观看| 天天碰天天摸天天搞| 人妻系列级片在线观看视频| 99在线视频精品观看高| a级黄片免费观看| 国产午夜在线播放视频| 国产极品气质外围av| 亚洲三级综合在线观看| 日韩一区二区在线播放观看| 亚洲一区二区偷拍女厕所| 91亚洲最新蜜桃在线| 亚洲熟女乱色一区二区三区视频| 精品人妻人人做人人爽| 亚洲欧美另类校园春色| 午夜精品久久秘?18免费观看| 美国伦理片午夜理论片| 天天在线播放日韩av| 亚洲国产电影的一区| 美女福利视频一区二区三区四区| 成人十欧美亚洲综合在线 | 九十九步都是爱最后一步是尊严| 人妻少妇精品二三区| 99福利一区二区视频| 中文字幕久久久国产| 精品国产无乱码一区二区三区 | 96在线观看免费播放| 青娱乐这里只有精品| 久久久久夜色国产精品电影| 人人妻人人爽人人摸| 亚洲熟妇丰满多毛xxxx网站| 夜夜操天天干夜夜操| 老司机免费视频福利0| 天堂av国产av伦理av| 国产av精品一区二区三区久久| 呻吟求饶的人妻中文字幕| 91 精品视频在线看| 欧美精品激情在线不卡| 中文字幕综合网91| 日韩女同与成人用品电影免费看| 欧美成人少妇人妻精品| 第一福利视频在线观看| 91九色91在线视频| 久久精品国产亚洲av清纯| 午夜美女福利视频在线| 国产精品蝌蚪自拍视频| tushy一区二区三区视频| 最新日韩中文字幕啪啪啪| 69精品互换人妻4p| 精品国模一区二区三区欧美| 99精品久久精品一区二区| 五十岁熟女高潮喷水| 日本亚洲午夜福利一区二区三区| 亚洲熟女乱色一区二区三区视频| 豆豆专区操逼性视频在线| 久久久久久久久久久久久国产| 青青在线免费手机播放视频| 大陆中文字幕视频在线| 一区二区三区四区视频精品免费| 青青草原在线播放日韩| 成熟了的熟妇毛茸茸| 首页欧美日韩中文字幕| 18岁禁一二三区免费体验| 快进来插我的逼嗯啊视频| 全国熟妇精品一区二区免费视频| 夫亡人妻被强干中文字幕| 欧美一区二区三区爽爽| 国产igao激情在线视频入口| 中文字幕人妻一区色偷偷久久| 一二三四区国产在线观看| 日本福利网站一区二区| 98热视频精品在线观看| 国产极品气质外围av| 久久久亚洲熟女一区二区| 大奶熟妇激情操逼逼| 女人的天堂 av在线| 亚洲经典av中文字幕| 91精品国产人妻麻豆| 久久视频 在线播放| 欧美成人红桃视频在线观看| 38av一区二区三区| 国产精品 亚洲欧美 自拍偷拍| 最近最新欧美日韩精品| 爱搞视频在线观看视频91| 国产精品性感美女视频| 网友自拍第一页99热| 欧美第一激情综合网欧美激情| 亚洲熟妇丰满多毛xxxx网站| 天堂网免费在线电影| 青青操91美女国产| 久久久久久久岛国免费观看| 18在线观看免费观看| 亚洲精品国产99999| 日本人妻熟妇丰满成熟HD系列| 91精品国产成人久久久久久| 午夜国产精品免费视频| 果冻麻豆一区二区三区| 男女爱爱好爽视频免费看| 69久久夜色精品国产69乱电影| 欧美成人久久久桃色aa| 少妇熟女天堂网av| 日本熟妇乱妇熟色视频| 核xp工厂精品久久亚洲| 97成人老师在线视频| 日本黄页在线观看视频| 亚洲第一中文字幕成人| av中文字幕国产精品| 可以直接看av网站| 久久无码高清免费视频| lutu玩弄人妻短视频| 欧美视频免费观看777| 2020年亚洲男人天堂网| 国产精品福利久久久久| 亚洲欧美国产人成在线| 欧美黑人1区2区3区| 91精品国产91久久久久久密臀| 91佛爷视频在线观看| 五月天男人的天堂中文字幕 | 亚洲精品1卡2卡3卡| 久草视频在线视频在线视频| 99久9在线视频播放| av激情四射五月婷婷| 国产天堂av不卡网| 中文字幕一区二区三区久久久| 可以免费观看日韩av| 亚洲日本欧美韩国另类综合| 麻豆出品视频在线观看| av人摸人人人澡人人超碰小说| 韩国一级片最火爆中文字幕| 老熟女xxxⅹhd老熟女性| 日本a级2020在线观看| 50熟妇一区二区三区| 一区二区三区国产在线成人av| 一区二区三区内射美女| 亚洲综合在线视频在线播放| 亚洲国产电影的一区| 99 re国产精品| 亚洲熟妇丰满多毛xxxx网站| 国产精品网站的黄色| 国产做A爱免费视频在线观看| 成年人免费黄色av| 伊人网在线欧美日韩在线| 日本高清有码在线视频| 68福利精品在线视频| 亚洲成人中文无码在线| 91精品国产欧美在线| 国产精品国产三级在线高清观看 | 精品国模一区二区三区欧美| 国产精品性感美女视频| 午夜久久人妻一级内射av网址 | 亚洲综合另类欧美久久| 日本高清在线观看不卡视频| 国内精品一区二区2021在线| alisontyler和黑人| 中文字幕国产一区在线视频| 国产资源在线观看二区| 青青青青午夜手机国产视频| 涩涩黄片在线免费观看| 夜夜操天天干夜夜操| 91精品国产91久久久久久密臀| 99国产精品久久99久久久| 中文字幕日韩人妻在线三区| 国产 亚洲 欧美 自拍| 神马不卡视频在线视频| 91美女在线观看视频| 天天插天天干天天狠| 久久人妻诱惑我视频| 成人av中文字幕在线看| 在线看的免费网站黄| 五月天男人的天堂中文字幕| 熟妇高潮久久久久久久| 日韩av水蜜桃一区二区三区| 人妻少妇视频系列视频在线| 久久99精品热在线观看| 在线看的免费网站黄| 欧美一区日韩二区三区四区| 开心五月综合激情婷婷| 操烂你的骚逼天天欧美| 97成人老师在线视频| 日本福利视频网站导航| 97超碰人人爽人人做| 久草视频在线视频在线视频| 91精品麻豆91夜夜骚| 91久久久久久最新网站| 精品不卡一区二区三区| 人妻色综合aaaaaa网| 农村大炕有肉大屁股熟妇| 91在线九色porny| 91porny九色视频偷拍| 日本黄页在线观看视频| 黄色网络中文字幕日本| av里面的动作是真进去吗| 亚洲韩精品一区二区三区| 美国十次了亚洲天堂网国产| 九一精品人妻一区二区三区| 美女网站福利在线观看| 亚洲综合另类欧美久久| 七色福利视频在线观看| xxoo福利视频导航| 新香蕉视频香蕉视频2| 午夜情色一区二区三区| 亚洲国产综合久久精品| 丝袜美女诱惑佐佐三上 | 男人电影天堂在线观看| 国产资源网站在线播放| 免费啪啪啪网站在线观看| 日本黄页在线观看视频| 得得爱在线视频观看| 一区二区三区四区视频精品免费| 中文字幕人妻一区色偷偷久久| 999国产精品视频免费看| 一二三四区国产在线观看| 加勒比不卡在线视频| 999国产精品视频免费看| 天天日夜夜操人人爽| 亚洲人人爽人人澡起碰av| 麻豆国产精品777777在| 玖辛奈18禁同人污本子| 欧美精品999不卡| 美女黄色啊啊啊啊视频| 夜夜人人干人人爱人人操| av在线播放观看h| 97人妻人人揉人人躁人人夜夜爽| 亚洲激情视频在线观看免费| 日本熟妇乱妇熟色视频| 亚洲av日韩久久网站| 91激情四射婷婷综合| 日本少妇人妻中文在线| 午夜精品一区二区三区不卡顿| 精产国品一二三77777| 日本少妇三级交换做爰做| 久久亚洲国产成人精品麻豆 | 人妻少妇精品二三区| 熟妇高潮久久久久久久| 蜜桃tv一区二区三区| 成熟了的熟妇毛茸茸| 亚洲va999天堂va| 人妻激情综合久久久久蜜桃| 国产 亚洲 欧美 自拍| 免费的啪啪视频软件| av资源中文字幕在线观看| 亚洲成人五月婷婷久久综合| 亚洲中文字幕最新地址| 熟妇人妻丰满久久久久久久| 国产美女主播av在线| caopeng97在线观看视频| 国色天香一二三期区别大象| 日本高清有码在线视频| 黄色av网址在线播放| 亚洲熟女人妻自拍在线视频 | 荣立三等功退休有什么待遇| 秋霞成人午夜鲁丝一区二区三区| 午夜国产免费视频亚洲| 亚洲午夜精品一级毛片app| 亚av一二三在线观看| 一区二区在线观看视频观看| 青青免费观看视频| 成人精品影视一区二区| 69视频在线精品国自产拍| 高清国产美女a一级毛片| 91性高湖久久久久久久久久| 人妻熟女 亚洲 一页二页| av资源中文字幕在线观看| 亚洲中文字幕在线视频观看二区| 亚洲男人的天堂最新网址| 日本不卡视频一二三区| 污网址在线观看视频| 女人扒开逼让男人操| 3344永久在线观看视频下载| 外国美女舔男人坤坤| 河北全程露脸对白自拍| 91九色人妻在线播放| 久久人妻诱惑我视频| 亚洲欧美另类丝袜另类自拍| 国产男人的天堂一区| 日本熟女0930视频| 美女av色播在线播放| 美女把逼扒开让男人桶| 欧美日韩精品aaa| 亚洲一区二区三区无码在线| 亚洲av手机免费在线| 网友自拍第一页99热| 天天干天天日天天弄| 欧美日韩高清片在线观看| 人妻免费视频黄片在线视频| 91精品麻豆91夜夜骚| 91精品资源在线观看| 亚洲成人激情在线综合| 日韩精品视频一区二区三区在线| 精久久久久久久久久久久| 内地精品毛片在线观看| 制服丝袜中文字幕熟女人妻| 人妻人妻在线视频网站| 亚洲成人自拍av在线| 美女张开腿给男人桶爽的软件| 琪琪日本福利伦理视频| 日韩一区二区在线播放观看| 免费在线观看亚洲福利| 操死你美女在线视频| 最近最新最好看的中文字幕| av毛片在线观看网址| 男女爱爱好爽视频免费看| 精品美女洗澡一区二区| av天堂新资源在线| 九九视频在线观看全部| 欧美操大黑鸡巴视频在线观看| 国产精品成人免费电影| 中文字幕 一区二区在线观看| 欧美vs亚洲vs日韩| 美女网站视频久久精品| 97超碰人人爽人人做| 青青在线视频看看| 国产91黑丝小视频在线观看 | 欧美插插插插插插| 川上优所有中文字幕在线| 老司机伊人99久久精品| 自拍偷拍 亚洲性图 欧美另类| 亚洲国产日韩a在线欧美| 国产黑色丝袜 在线日韩欧美| 久久午夜免费鲁丝片| 国产激情一区二区视频| 中文字幕精品人妻久久久久| 77亚洲视频在线观看| 抽插小穴啊啊啊视频| 日本高清久久人人爽| 五月激情婷婷四射基地| 港台美女明星av天堂| 国产又粗又长又大视频| 亚洲综合在线视频在线播放| 亚洲天堂av最新在线| 午夜野花视频在线观看| 免费在线观看视频啪啪| 亚洲高清免费在线观看视频| 久久99久久99久久97的人| 亚洲a区在线免费观看| 欧美视频免费观看777| 亚洲男人天堂最新网址大全 | www国产亚洲精品久久久| 99久久国产精品免费消防器材| www,日韩av,com| 久久久西西gogo日本美女人体 | 欧美大鸡吧男操女啊啊啊视频| 五月天色婷婷狠狠爱| 欧美极品少妇高潮喷水| 国产精品免费看一区二区三区| 久久国产半精品99精品国产| 人妻被强av系列一区二区| 天天日天天玩天天摸| 最新国产午夜激情视频| 亚洲字幕一区二区夜色av| 欧美猛少妇色ⅹⅹⅹⅹⅹ猛叫| 九九热视频1这里只有精品| 欧美最新一区二区三区| 国产人妻777人伦精品hd超碰| 啊不行啊操逼好爽大鸡吧视频| av资源中文字幕在线观看| 美女网站视频久久精品| 午夜精品秘一区二区三区| 老牛影视在线一区二区三区| 538欧美在线观看一区二区三区| 国产精品美女免费视频观看| 插鸡视频免费网站在线播放| 狂操鸡巴小骚逼视频免费观看| 熟妇高潮久久久久久久| 四虎精品久久免费最新| 亚洲永远av在线播放| 亚洲精品1卡2卡3卡| av无限看熟女人妻另类av| 91超碰国产在线观看| 成人资源中文在线观看| 狠狠操狠狠操狠狠插| 午夜精品一区二区三区不卡顿| 午夜92福利1000| 天天日天天干天天日天天干天天| 亚洲成a人77777| 欧美熟女xx00视频| 北野中文字幕一区二区| 亚洲成人av在线一区二区| 久久国产半精品99精品国产| 亚洲午夜精品视频节目| 一区二区三区观看在线| 日本一区二区高清av中文| 亚洲欧美综合另类最新| 欧美成人久久久桃色aa| alisontyler和黑人| 69精品互换人妻4p| xxxx69在线观看视频| 五十岁熟妇高潮喷水| 黄在线看片免费人成视频| 最新国产精品久久精品app| 最新日韩中文字幕免费在线观看| 在线观看网站伊人网| 69精品互换人妻4p| 91进入蜜桃臀在线播放| 国际日韩日韩日韩日韩日韩 | 丝袜美腿日韩av一区| xxnxx国产美女| 中文字幕福利视频在线一区| 女人的天堂av在线网| 欧美亚洲另类精品第一页| 成人人妻h在线观看| 欧洲成熟女人色惰片| 久久久国产精品免费视频网| 高潮喷水一区二区三区| 亚洲一区视频中文字幕在线播放 | 亚洲欧美日韩电影一区| 少妇被中出一区二区| 在线成人教育平台排名| 亚洲综合第一区二区| 欧美成人屋影院在线视频观看| av福利免费体验观看| 97香蕉久久国产超碰| 亚洲图片另类综合小说| 中文字幕 首页 人妻| 欧美日韩福利视频网| 欧美人与动欧交视频| 亚洲无人区乱码中文字幕一区| 中文字幕综合网91| 强乱人妻中文字幕日本| 插鸡视频免费网站在线播放| 久久久久高潮白浆久久| 在线免费视频999| 最新日韩av电影在线播放 | 一区二区三区高清视频3| 不卡高清一区二区三区| 日韩欧美一区二区三区免费看| 午夜宅男电影av网站| 蜜乳视频一区二区三区| 欧美啪啪一区二区三区| 抽插小穴啊啊啊视频| 91大神福利视频网| 中文字幕欧美人妻在线.| 大尺度久久久久久久| 在线免费观看视频18| 免费看一级高潮喷水片| 人妻激情偷乱一区二区三区av| 丝袜美女诱惑佐佐三上| 亚洲天堂色综合久久| 亚洲自拍偷拍一区二区中文字幕| 欧美成人性生活视频播放| 5d蜜桃臀女无痕裸感| 日韩国产欧美久久一区| 一区二区三区午夜福利在线| 国产成人情侣激情视频| 中文字幕在线免费观看人妻| 麻豆白洁少妇在线播放| 欧美一级特黄大片在线| 国产精品美女免费视频观看 | avgo成人短视频| 人妻熟女 亚洲 一页二页| 欧美日韩一区二区三区成人影院| 国产av在线免费视频| 亚洲精品综合欧美精品综合| 成人十欧美亚洲综合在线| 亚洲乱码av一区二区蜜桃av| 欧美在线观看一区二区不卡| 黄版视频在线免费观看| 国产精美视频精品视频精品| av成人三级高清日韩| 国产成人情侣av在线| 嗯~嗯~啊啊啊~高潮了软件| 中文字幕福利视频第四页| 国产精品内射婷婷一级| 亚洲国内精品久久久久久久| 午夜精品小视频在线播放| 亚洲中文字幕无线乱码人妻精品| 在线免费观看a视频免费| 99国产精品国产精品毛片19| 国产最新av在线免费观看| 欧美黄色性视频网站| 1区3区4区产品乱入视频| 天天早上头和脸出汗是怎么办| 国产精品免费看一区二区三区| 国产精品成人免费电影| av在线中文字幕在线| 在线观看视频免费一区二区三区| 狂操鸡巴小骚逼视频免费观看| av激情四射五月婷婷| 东京热男人的天堂视频| 一区二区在线观看视频观看| av成人三级高清日韩| 日韩av水蜜桃一区二区三区| 青青在线视频看看| 最近日韩免费在线观看| 最新国产午夜激情视频| 最近在线中文字幕免费| 有码一区二区三区四区五区| 黄色网络中文字幕日本| 国产av在线免费视频| 久草视频在线视频在线视频| a级片特黄免费看| 老司机在线视频福利观看| 日韩人妻中文字幕区| 青娱乐这里只有精品| 99国产精品国产精品毛片19| 免费高清av一区二区| 天堂网成人av电影| 中文字幕熟女人妻丝袜丝在线| 成人超碰一区二区三区| 国产视频成人一区二区| 亚洲一区二区偷拍女厕所| 国产午夜在线播放视频| 美国十次了亚洲天堂网国产| 午夜福利午夜福利影院| 亚洲天堂色综合久久| 精品一区二区三区喷水内射高潮| 一区二区三区四区 在线播放 | 91在线九色porny| 国产在线小视频一区二区| 日本韩国福利在线播放| 亚洲女人自熨在线视频| 久久精品久久久久观看99水蜜桃| 天天日天天亲天天操| av成人三级高清日韩| 人妻被强av系列一区二区| 亚洲综合一区二区三区四区| 亚州av嫩草av极品在线观看| 99久久国产精品免费热| 人妻超清中文字幕在线乱码| 在线免费观看视频18| 91精品夜夜夜一区二区蜜桃| 自拍偷自拍亚洲精品10p| 97超碰人人爽人人做| 中文字幕在线观看av观看| 国产91精品福利系列| 人妻中文字幕亚洲在线| 999国产精品视频免费看| 国产美女视频带a∨黄色片| 亚洲一区视频中文字幕在线播放| 漂亮人妻口爆久久精品| 国产乱码有码一区二区三区| 三级欧美日韩一区二区三区| 91超碰九色porny| 亚洲欧美日韩电影一区| 日本不卡视频一二三区| 天天插天天操天天射天天干| 3344永久在线观看视频下载| 一区二区三区四区影片| 加勒比东京热绿帽人妻多人操| 青青青在线视频免费播放| 啊不行啊操逼好爽大鸡吧视频| 天堂av国产av伦理av| 深夜福利免费观看在线看| 久久久久高潮白浆久久| 亚洲欧美韩国日本一区二区| 青青草一个释放的网站| 久久精品久久久久观看99水蜜桃| 亚洲国产精品久久久久久无码| 亚洲熟女少妇中文字幕系列| 中国精品人妻一区二区| 国产不卡免费在线观看| 国产精品午夜无码AV体验区| 91九色人妻在线播放| 九九六视频,这里只有精品| 国产精品成人免费电影| 天天插天天干天天狠| 亚洲无码专区中文字幕专区| 午夜野花视频在线观看| 91九色91在线视频| 天天干天天色综合久久| 人妻系列中文字幕大乳丰满人妻| 中文字幕观看中文字幕免费| 天天摸天天干夜夜操| 亚洲国产精品一区二区第二页| 中日韩又粗又硬又大精品| 久久久久久久精品乱码| 亚州av嫩草av极品在线观看| 亚洲精品乱码久久久久app| 2018中文字字幕人妻| 人妻激情综合久久久久蜜桃| 欧美色视频网址大全| 亚洲人人爽人人澡起碰av| 男人的天堂aⅴ在线| 国产精品剧情av在线播放| 国产igao激情在线视频入口| av福利免费体验观看| 免费中文字幕a级激情| 欧美第一激情综合网欧美激情| 自拍偷拍亚洲综合第一页| 亚洲另类欧美综合久久| 99999久久久精品| 日本少妇丰满大bbb的小乳沟| 东北老女人熟女啪啪视频| 91人妻人人做人人爽高清| 大陆中文字幕视频在线| 天天干天天色综合久久| 99久久国产精品免费热| 亚洲av三级电影在线观看| 午夜情色一区二区三区| 狠狠操狠狠操狠狠插| www一区二区91| 精品免费一区二区三区四区视频 | 久久人妻诱惑我视频| 午夜国产免费视频亚洲| 亚洲精品国产99999| 精品不卡一区二区三区| 日韩精品视频一区二区三区在线| 国产精品网站亚洲发布| 神马午夜久久电影网| 国产精品视频网站污污污| 女人的天堂av在线网| 开心激情五月天作爱片| 欧美日韩一区二区三区成人影院| 午夜野花视频在线观看| 国产在线小视频一区二区| 青青青在线观看国产| 偷拍欧美日韩另类图片| 美女福利视频一区二区三区四区 | 91精品夜夜夜一区二区| 黑人和日本人av一区二区| 青青青在线视频观看97| 日本香港韩国三级黄色| 一区二区三区国产精华液区别大吗| 午夜一区二区三区视频在线观看| 加勒比东京热绿帽人妻多人操| 精品日本少妇久久久|