版块导航
正在加载中...
客户端APP下载
论文辅导
登录
注册
帖子
帖子
用户
本版
应《网络安全法》要求,自2017年10月1日起,未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用,请尽快对帐号进行手机号验证,感谢您的理解与支持!
24小时热门版块排行榜
>
论坛更新日志
(2698)
>
虫友互识
(170)
>
找工作
(75)
>
考博
(47)
>
导师招生
(39)
>
基金申请
(39)
>
硕博家园
(39)
>
职场人生
(38)
>
文献求助
(23)
>
招聘信息布告栏
(22)
>
休闲灌水
(15)
>
考研
(10)
>
育儿交流
(9)
>
论文投稿
(8)
>
专利求助
(7)
>
标准与专利
(7)
小木虫论坛-学术科研互动平台
»
计算模拟区
»
分子模拟
»
Monte Carlo
»
【转帖】2维无外场Ising模型的Monte Carlo模拟的C++实现
8
1/1
返回列表
查看: 2418 | 回复: 7
只看楼主
@他人
存档
新回复提醒
(忽略)
收藏
在APP中查看
【奖励】
本帖被评价4次,作者zyj8119增加金币
3
个
zyj8119
木虫
(著名写手)
模拟EPI: 10
应助: 65
(初中生)
贵宾: 0.003
金币: 1035.1
帖子: 2936
在线: 1329小时
虫号: 664177
[
资源
]
【转帖】2维无外场Ising模型的Monte Carlo模拟的C++实现
CODE:
// Copyright(c)2010 Ma Chengzhang
#include
#include
#include
//tr1随机数调用
#include
//设置随机数种子用
#include
using namespace std; //cout,endl等简单语句必须
using namespace std::tr1; //tr1随机数调用
inline unsigned __int64 ReadTSC()//为取到精确的随机数种子,嵌入汇编程序取系统时钟
{
__asm _emit 0x0F
__asm _emit 0x31
}
double totalE(vector
>& SpinGrid,int GridL,double JBound)
{
int nn_count = 0;//最近邻计数
// 每个格点只考虑右键和上键
for(int i=0; i
{
for(int j=0; j
{
if(SpinGrid[i][j]==SpinGrid[i+1][j]) // 考虑右键
nn_count++;
else
nn_count--;
if(SpinGrid[i][j]==SpinGrid[i][j+1]) // 考虑上键
nn_count++;
else
nn_count--;
}
}
//周期性边界情况
for(int i=0; i
{
if(SpinGrid[i][GridL-1]==SpinGrid[i+1][GridL-1]) // 右键
nn_count++;
else
nn_count--;
if(SpinGrid[i][GridL-1]==SpinGrid[i][0]) // 上键
nn_count++;
else
nn_count--;
}
for(int j=0; j
{
if(SpinGrid[GridL-1][j]==SpinGrid[0][j]) // 右键
nn_count++;
else
nn_count--;
if(SpinGrid[GridL-1][j]==SpinGrid[GridL-1][j+1]) // 上键
nn_count++;
else
nn_count--;
}
//右上顶点
if(SpinGrid[GridL-1][GridL-1]==SpinGrid[0][GridL-1]) // 右键
nn_count++;
else
nn_count--;
if(SpinGrid[GridL-1][GridL-1]==SpinGrid[GridL-1][0]) // 上键
nn_count++;
else
nn_count--;
return -JBound*nn_count;
}
void main()
{
//物理常数定义
double kB=1.0; //Boltzmann常数
double JBound=0.5; //交换作用系数
int GridL=20; //网格大小,长宽一样,即 GridL x GridL 网格
mt19937 RandEngine; //Mersenne Twister算法随机数生成器,核心引擎类
unsigned long RandSeed = (unsigned long)ReadTSC(); //用 CPU 时钟作为随机数种子
RandEngine.seed(RandSeed);
uniform_int
RandInt(0, GridL-1);//产生0到GridL-1之间的随机整数,
// 含0和GridL-1,即[0,GridL-1]而不是(0,GridL-1)
uniform_real
RandFloat(0,1); //产生0到1之间的随机数
double MaxTemp=5.00; //要模拟的最高温度,温度以J/kB为单位,kB是波尔兹曼常数
double MinTemp=0.1; //要模拟的最低温度,温度以J/kB为单位,kB是波尔兹曼常数
unsigned int BeforeEquilSteps=10000; //在达到平衡前要演化的步数
unsigned int TotalSteps=BeforeEquilSteps+10000; //每温度要演化的总步数,包含前面的BeforeEquilSteps
int DiscardedSteps=50; //达到平衡态后,为实现无偏统计抽样,每DiscardedSteps演化步数才统计一次
double TempInterval=0.1; //扫描温度间隔,温度以J/kB为单位,kB是波尔兹曼常数
double beta; //当前beta=1/(kB*T)
vector< vector
> SpinGrid(GridL, vector
(GridL,-1)); //网格定义,格点初始值都取-1
//统计量定义
double TotalEnergy=0; //网格总能量
double persiteEnergy=0; //每格点能量
double persiteAveEnergy=0; //每格点平均能量
ofstream fileout("Ising2D.dat");
//Monte Carlo 部分
for(double CurrentTemp=MinTemp;CurrentTemp
{
int count=0; //此变量记录达到平衡后演化步数,每当count达到DiscardedSteps时归零
int Outcount=0; //记录累计统计次数
int left,right,bottom,top; //辅助变量
beta=1/(kB*CurrentTemp);
for(unsigned int mcStep=0;mcStep
{
//随机选择格点反转自旋,随机数算法取自:C++ tr1 Mersenne Twister算法
int i=RandInt(RandEngine); //第i列
int j=RandInt(RandEngine); //第j行
//采用周期性边界条件
if(i!=0)
left = SpinGrid[i-1][j];
else
left = SpinGrid[GridL-1][j];
if(i!=GridL-1)
right = SpinGrid[i+1][j];
else
right = SpinGrid[0][j];
if(j!=0)
bottom = SpinGrid[i][j-1];
else
bottom = SpinGrid[i][GridL-1];
if(j!=GridL-1)
top = SpinGrid[i][j+1];
else
top = SpinGrid[i][0];
int spinSum=top+bottom+left+right;
//double oldEnergy=-JBound*SpinGrid[i][j]*spinSum;
//double newEnergy=-oldEnergy;
double EDiff=2*JBound*SpinGrid[i][j]*spinSum;//EDiff=newEnergy-oldEnergy
//cout<<"EDiff="<
//判断是否反转自旋
if( EDiff<0)
SpinGrid[i][j]=-SpinGrid[i][j];
else
{
if(RandFloat(RandEngine)
SpinGrid[i][j]=-SpinGrid[i][j];
}
//如果mcStep达到了指定步数,就认为达到了平衡态,可以进行状态分析了
//但是只是每隔DiscardedSteps步数才分析一次
if(mcStep>BeforeEquilSteps)
{
count=count+1;
if(count==DiscardedSteps)
{
count=0; //count归零,重新开始计数
TotalEnergy=totalE(SpinGrid,GridL,JBound);
persiteEnergy=TotalEnergy/(GridL*GridL);
persiteAveEnergy=(persiteAveEnergy*Outcount+persiteEnergy)/(Outcount+1);
Outcount++;
}
}
}
//输出数据到数据文件
fileout<< CurrentTemp <<" "<
//统计量归零
persiteAveEnergy=0;
Outcount=0;
}
//Monte Carlo 部分结束
fileout.close();
////输出测试
////测试随机整数
//for(int i = 0; i < 50; ++i) //产生均匀分布的在0到19之间的50个整数随机数
//{
// std::cout << RandInt(RandEngine) << endl;
//}
//cout<<"*************"<
//for(int i = 0; i < 50; ++i) //产生均匀分布的在0到1之间的50个随机数
//{
// std::cout << RandFloat(RandEngine) << endl;
//}
////测试随机整数毕
//cout<
//cout<<"TotalSteps = "<
//cout<<"RandSeed = "<
///////////////////////////////////////////////////////////////////////
cout<
char c;
cin>>c;
}
下面是用Mathematica处理数据的结果:
[
Last edited by zyj8119 on 2010-9-10 at 17:51
]
回复此楼
» 收录本帖的淘贴专辑推荐
资源收集
graphene
学习
量子蒙特卡洛模拟
» 猜你喜欢
硕博巨婴,也许才刚刚开始
已经有69人回复
杭电、天津科技、青农和宁波工程学院如何选?
已经有20人回复
34岁读博士晚吗
已经有26人回复
初始合伙人来啦!(生物试剂耗材标准品)
已经有6人回复
江西双非一本和四川双一流高校如何选择?
已经有23人回复
投稿被一个审稿人恶意评审了怎么样?
已经有7人回复
面青地会评时间
已经有5人回复
关于2023的收入
已经有34人回复
累并快乐着
已经有13人回复
每次骚扰女学生的都是院系领导,而不是普通教师,小编们要注意措辞正确
已经有9人回复
» 本主题相关价值贴推荐,对您同样有帮助:
泵站进水流道出口断面轴向流速分布不均匀度的求解方法
已经有4人回复
matlab解微分方程
已经有10人回复
模拟结果实验结果不一致
已经有11人回复
fluent中的体热源加载问题(实际上是热流密度,不是热生成率)
已经有19人回复
dynaform 软件模拟问题,三通管成形的冲头怎么设置
已经有5人回复
纤维素酸解的相关资料
已经有4人回复
C++编程,关于循环结构的,大家看看我这程序哪里错了?
已经有21人回复
Ansys Fluent的计算能力有多大(最多能计算多少网格?)
已经有18人回复
用DPM模型时,solution里面会不会有 discrece model方程
已经有4人回复
是什么原因导致纤维素的耐盐性不好?
已经有6人回复
Monte Carlo 模拟时图是用什么软件做的
已经有10人回复
【求助】羟丙基纤维素和聚乳酸的共同溶剂
已经有6人回复
【求助】用monte carlo法做分子模拟,用什么编程语言实现比较好?
已经有25人回复
【原创】抛土引玉, Monte Carlo 如何自己写程序: NVT_LJ
已经有80人回复
【讨论】大家对未来GPU在monte carlo中的应用持有什么样的态度(CUDA)?
已经有15人回复
【求助】有用Monte Carlo 编过有关Fluid中粒子运动的程序
已经有6人回复
【求助】求Monte Carlo 关于扩散系数计算的源代码
已经有5人回复
【在线答疑】经典粒子体系的Monte Carlo 模拟之基础篇
已经有73人回复
1楼
2010-09-10 17:49:48
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
yahoohoo
铁杆木虫
(著名写手)
模拟EPI: 10
应助: 0
(幼儿园)
贵宾: 1.55
金币: 7632.4
帖子: 1176
在线: 167小时
虫号: 74894
★ ★ ★ ★ ★ ★ ★
zh1987hs(金币+5):谢谢yahoohoo的精彩回复 2010-09-10 22:54:51
lei0736(金币+2):谢谢 2010-09-11 08:44:08
引用回帖:
Originally posted by
zyj8119
at 2010-09-10 17:49:48:
[code]// Copyright(c)2010 Ma Chengzhang <
Machengzhang@126.com
>
#include
#include
#include
//tr1随机数调用
#include
//设置随机 ...
这是你自己写的程序吗?很不错啊,:-)。
这是个巨正则系综下的Monte Carlo模拟,因为总的自旋 $ S = \sum_i {s_i} $ 是变化的。
随机数的种子是无所谓精确与否的,代码中的做法只是为了使得不同的模拟使用不同的随机数序列。
建议将程序模块话,比如 Ising MC的一般模块为:
void INIT(); // Set initial configuration
double totalEnergy(); // Calc total energy
double SpinEnergy(int i, int j); // Energy of a given spin
void flip(); // flip a randomly chosen spin
主程序部分的周期性边界条件处理可以更快:
int iRight = i + 1, iLeft = i - 1, jUp = j + 1, jDown = j - 1;
if (iRight = BOXLX) iRight = 0;
if (iLeft = -1) iLeft = BOXLXM1; // BOXLXM1 = BOXLX - 1;
if (jUp = BOXLY) jUp = 0;
if (jDown = -1) jDown = BOXLYM1; // BOXLYM1 = BOXLY - 1;
这样做的好处是将if else 分支结构转变成单独的if, 既形式简单又速度更快
能量判断可以简单写成如下代码:
double dE = Enew - Eold;
if (dE <= 0 || rng.gen_open0_open1() < exp(-BETA * dE)) { // 这样的好处是如果dE <= 0 成立那么后面的部分就不用计算
// ACCEPTED
}
else {
// REJECTED
}
关于能量与温度的关系,建议能量用 $ E/NJ $, 这样的话在低温区,所有的自旋方向一致,那么构象能为 $ E_{min} = -2 N |J| $。所以最低能量应该为 $ E / NJ = -2 $。而在高温区,自旋向上与向下的概率相等,因为每个自旋近邻的四个自旋和为0,因为总能量也为0。根据这两点,模拟得到的结果是有问题的。
[
Last edited by yahoohoo on 2010-9-10 at 21:13
]
赞
一下
(2人)
回复此楼
2楼
2010-09-10 21:05:59
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
zyj8119
木虫
(著名写手)
模拟EPI: 10
应助: 65
(初中生)
贵宾: 0.003
金币: 1035.1
帖子: 2936
在线: 1329小时
虫号: 664177
引用回帖:
Originally posted by
yahoohoo
at 2010-09-10 21:05:59:
这是你自己写的程序吗?很不错啊,:-)。
这是个巨正则系综下的Monte Carlo模拟,因为总的自旋 $ S = \sum_i {s_i} $ 是变化的。
随机数的种子是无所谓精确与否的,代码中的做法只是为了使得不同的模拟使 ...
谢谢建议!!
回复此楼
3楼
2010-09-10 22:56:11
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
zyj8119
木虫
(著名写手)
模拟EPI: 10
应助: 65
(初中生)
贵宾: 0.003
金币: 1035.1
帖子: 2936
在线: 1329小时
虫号: 664177
引用回帖:
Originally posted by
yahoohoo
at 2010-09-10 21:05:59:
这是你自己写的程序吗?很不错啊,:-)。
这是个巨正则系综下的Monte Carlo模拟,因为总的自旋 $ S = \sum_i {s_i} $ 是变化的。
随机数的种子是无所谓精确与否的,代码中的做法只是为了使得不同的模拟使 ...
不是我写的,我是转帖的。
赞
一下
回复此楼
4楼
2010-09-11 12:02:55
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
咔咔咔咔
新虫
(初入文坛)
应助: 0
(幼儿园)
金币: 11.6
帖子: 35
在线: 26.3小时
虫号: 1348572
★★★★★ 五星级,优秀推荐
用vc6.0可以不
回复此楼
5楼
2012-04-26 19:13:43
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
yaneehan2010
银虫
(正式写手)
应助: 8
(幼儿园)
金币: 548.9
帖子: 420
在线: 106.3小时
虫号: 1044536
★★★★★ 五星级,优秀推荐
你好,我是刚接触分子模拟的菜鸟,做的是气体分子扩散系数和分子间碰撞情况。我该怎么着手呢?不知道有没有类似的程序可以下载呢?在哪里可以下载的到呀,请指条明路,万分感谢、、
赞
一下
回复此楼
6楼
2012-05-23 00:43:45
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
马pi精
新虫
(小有名气)
应助: 0
(幼儿园)
金币: 32.4
帖子: 69
在线: 17.9小时
虫号: 1417222
★★★ 三星级,支持鼓励
楼主试过画出T在0~1内的比热曲线吗? 求交流
赞
一下
回复此楼
7楼
2012-08-17 23:18:37
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
quantumage
金虫
(小有名气)
应助: 0
(幼儿园)
金币: 1380.8
帖子: 292
在线: 136.1小时
虫号: 1855704
★★★★★ 五星级,优秀推荐
给介绍几本量子monet carlo的学习书籍呗!我也想自己算一下Ising模型,和XY模型!
赞
一下
回复此楼
8楼
2013-01-09 00:15:45
已阅
回复此楼
关注TA
给TA发消息
送TA红花
TA的回帖
相关版块跳转
第一性原理
量子化学
计算模拟
分子模拟
仿真模拟
程序语言
我要订阅楼主
zyj8119
的主题更新
8
1/1
返回列表
☆ 无星级
★ 一星级
★★★ 三星级
★★★★★ 五星级
如果回帖内容含有宣传信息,请如实选中。否则帐号将被全论坛禁言
普通表情
龙
兔
虎
猫
高级回复
(可上传附件)
SCI英文润色翻译-EnPapers美团队-先修改后付款
百度网盘
|
360云盘
|
千易网盘
|
华为网盘
在新窗口页面中打开自己喜欢的网盘网站,将文件上传后,然后将下载链接复制到帖子内容中就可以了。
最具人气热帖推荐
[查看全部]
作者
回/看
最后发表
[
找工作
]
初始合伙人来啦!(生物试剂耗材标准品)
+3
欢快的小科研人
2024-06-15
6/300
2024-06-16 07:13
by
小李小李小李
[
论文投稿
]
投稿被一个审稿人恶意评审了怎么样?
+5
1chen
2024-06-14
7/350
2024-06-15 23:15
by
xy66xy
[
教师之家
]
关于2023的收入
+33
小龙虾2008
2024-06-10
34/1700
2024-06-15 23:01
by
zeolitess
[
文学芳草园
]
累并快乐着
+13
MYHLD521
2024-06-14
13/650
2024-06-15 22:59
by
zeolitess
[
基金申请
]
关于博后基金的bug问题
+6
lxr1991
2024-06-14
9/450
2024-06-15 21:17
by
since—2010
[
基金申请
]
希望今年自己国自然面上项目和老婆青年项目能中!
+5
恐龙爸爸
2024-06-14
5/250
2024-06-15 20:36
by
wanghuawei
[
教师之家
]
请问事业编制和年薪制冲突吗?
+6
ZHONGWU_U
2024-06-14
6/300
2024-06-15 20:16
by
Ermito
[
考博
]
希望能25博士入学,可提前一年做科研助理
+4
干饭版小太阳
2024-06-09
8/400
2024-06-15 18:36
by
独孤老狗
[
博后之家
]
山东大学(青岛)“天然药物生物智造”课题组 招聘“博士后”(年薪20.4-55.6万元)
+3
第二种态度
2024-06-11
5/250
2024-06-15 17:56
by
小懂事k
[
论文投稿
]
投了一篇4区的SCI,审稿人一个拒稿,一个小修,编辑给了大修。
+9
安稳22123
2024-06-13
10/500
2024-06-14 23:45
by
jurkat.1640
[
考研
]
物理化学一对一辅导
+3
林大diao
2024-06-12
5/250
2024-06-14 20:57
by
林大diao
[
硕博家园
]
关于硕博连读的一些疑问?
+4
Lwenter
2024-06-14
4/200
2024-06-14 14:32
by
ou0551
[
有机交流
]
ππ堆积会发生在有机溶剂中吗
5
+3
zibuyu0420
2024-06-13
4/200
2024-06-14 14:17
by
小肉干
[
论文投稿
]
最近写了一篇控制优化领域的文章,可以投哪里啊?有没有水一些的期刊推荐
+7
香瓜木香
2024-06-12
13/650
2024-06-14 07:05
by
香瓜木香
[
论文投稿
]
投稿后发现其他作者的邮箱填错了该怎么办呀
10
+4
在飞的猪
2024-06-13
6/300
2024-06-14 04:45
by
小虫子咔咔
[
基金申请
]
博士后面上项目状态还是专家评审吗
10
+9
Thatcheremu
2024-06-13
55/2750
2024-06-13 21:23
by
乌合麒麟
[
硕博家园
]
科研求助
+5
杲www
2024-06-12
6/300
2024-06-13 16:16
by
姓李名明
[
硕博家园
]
机械研究生如何拿到年薪40+w
+13
阿巴阿巴哦哦
2024-06-11
15/750
2024-06-13 15:40
by
113745685
[
论文投稿
]
with editor日期变更
+3
慎独的小花卷
2024-06-12
8/400
2024-06-13 11:00
by
慎独的小花卷
[
找工作
]
成都产品质量检测研究院
200
+3
鲸鱼663
2024-06-11
8/400
2024-06-13 08:10
by
加纳居士
信息提示
关闭
请填处理意见
关闭
确定