--- /dev/null
+ 0.02
+ -2.09251
+ -1.67129
+ -2.43220
+ -2.44579
+ -2.30991
+ -1.65770
+ -3.05723
+ -1.30442
+ -0.42122
+ 0.00000
+ -0.35328
+ 0.05435
+ 0.29893
+ 0.81526
+ 0.86961
+ 1.04625
+ -0.17664
+ 1.37236
+ 1.34518
+ -0.97831
+ -1.67129
+ -2.43220
+ -0.35328
+ -0.35328
.0572167103 -.0468608825 .0151048455 .0084963678
.0278930397 .0076922911 .1033536738 -.0098256036
.0611385674 .0448303346 .0861379100 .0861379100
+
+ 5.6053537261 6.2200032154 6.2159898496 6.4386968963
+ 6.2098101231 5.9596374798 5.4310222662 4.8790085257
+ 5.2638423822 5.4961954124 4.2928124162 4.3582250540
+ 4.2667755921 3.6813845958 3.5405429647 3.7026985912
+ 4.7177596156 3.3811992227 3.7936085458 4.5157490585
+ 6.2200032154 6.2159898496 5.2638423822 5.2638423822
+
+ 6.6296746680 6.6717260508 6.8168876945 6.8378960152
+ 6.3512039986 6.1425192423 5.5200206406 5.4858936024
+ 4.9382573262 4.2795129166 4.0534857018 4.2086428525
+ 3.5388086279 3.4209730327 2.7995744193 4.8164268553
+ 3.8790755071 3.5878262177 4.6457432555 6.6296746680
+ 6.6717260508 5.4858936024 5.4858936024
+
+ 6.6424306340 6.9715947250 6.9241225787 6.5291325353
+ 6.1394821460 5.4415971840 5.2914044780 4.7881880526
+ 4.1302408718 3.9405275117 3.9022731880 3.6078928145
+ 3.0809713206 3.0182595342 4.7935507272 3.8711928134
+ 3.6826168780 4.4872661030 6.6717260508 6.6424306340
+ 5.2914044780 5.2914044780
+
+ 7.3271707934 7.2970176793 7.0494948235 6.2565687662
+ 5.8133715113 5.8845382769 5.3463553241 4.7146928716
+ 4.4864793359 4.3092352779 3.6928361281 4.0100566871
+ 3.6584608489 4.5901676001 4.4681844010 4.6334334282
+ 4.9023945082 6.8168876945 6.9715947250 5.8845382769
+ 5.8845382769
+
+ 7.0106342493 6.8341607204 6.2047876082 5.5802122247
+ 5.6716968726 5.2011344932 4.1966522349 4.2383931955
+ 3.8954073453 3.5679086118 3.1136348127 2.7622561241
+ 4.5092053374 4.0238412606 3.9799853578 4.8087180377
+ 6.8378960152 6.9241225787 5.6716968726 5.6716968726
+
+ 6.5353606708 5.8844482033 5.1100976261 5.5515360893
+ 4.9587559238 4.2030790774 4.0351600709 3.6506289806
+ 3.4640416936 3.1611022087 2.7976192835 4.0422832763
+ 3.4415952740 3.8620694983 4.5934246009 6.3512039986
+ 6.5291325353 5.5515360893 5.5515360893
+
+ 5.2828811325 4.8298466315 4.7575327777 4.7249779042
+ 3.5656384285 3.5613366567 3.9273045564 3.7664779118
+ 3.5861699761 3.5523191627 4.6508672463 4.2424308260
+ 4.0987727783 4.5079669929 6.1425192423 6.1394821460
+ 4.7575327777 4.7575327777
+
+ 4.2222645754 4.1627149417 3.9786372270 3.1623319201
+ 2.9379696734 3.1365985188 2.8025851476 2.6576343487
+ 2.8733395837 3.8785891655 3.4443068252 3.4405132673
+ 4.0915305763 5.5200206406 5.4415971840 4.1627149417
+ 4.1627149417
+
+ 4.3860654351 3.6520014167 3.0264937881 2.5680579493
+ 2.5001857370 2.1393094782 1.4707077389 1.7373181916
+ 3.1045342752 2.0152329038 1.8196209136 3.6395852546
+ 5.4858936024 5.2914044780 4.3860654351 4.3860654351
+
+ 2.5016557935 2.7879606493 2.1424756816 1.0741544747
+ 1.3236438850 -.2712356678 .7613312181 2.8627348285
+ 1.5002127649 -.2392577636 3.5381128985 4.9382573262
+ 4.7881880526 3.6520014167 3.6520014167
+
+ 2.2903847447 2.2466771550 1.5053607856 1.4031797188
+ .9676209687 1.3075541842 2.8595118604 1.9833053538
+ -.0085773647 2.9557553666 4.2795129166 4.1302408718
+ 3.0264937881 3.0264937881
+
+ 1.2800038243 .7689806100 .9119606366 -.0744371793
+ .3669965427 2.4288563200 1.4975482174 -.2134101897
+ 2.9415021479 4.0534857018 3.9405275117 2.5680579493
+ 2.5680579493
+
+ -.6792428859 .4532383239 -.7590387660 -.3617034846
+ 1.6803275058 .6775209988 -.5354837468 2.6208591363
+ 4.2086428525 3.9022731880 2.5001857370 2.5001857370
+
+ .4136122500 -.1637152873 .2212564728 2.1287126020
+ .3896388275 -.2354089650 2.3255326077 3.5388086279
+ 3.6078928145 2.1393094782 2.1393094782
+
+ -3.5193531074 -2.0157307929 1.8355926160 2.7489174124
+ 2.0098391558 1.7975718667 3.4209730327 3.0809713206
+ 1.4707077389 1.4707077389
+
+ -1.3532688232 2.5504275249 2.8202873790 1.6253609671
+ 1.8625091247 2.7995744193 3.0182595342 1.7373181916
+ 1.7373181916
+
+ 3.7292778697 2.2944436481 -.0303565255 3.1115776177
+ 4.8164268553 4.7935507272 3.1045342752 3.1045342752
+
+ -.0827362961 -1.6043113182 2.4439837435 3.8790755071
+ 3.8711928134 2.0152329038 2.0152329038
+
+ -4.3897783280 2.3664634533 3.5878262177 3.6826168780
+ 1.8196209136 1.8196209136
+
+ 4.1927969260 4.6457432555 4.4872661030 3.6395852546
+ 3.6395852546
+
+ 6.6296746680 6.6717260508 5.4858936024 5.4858936024
+
+ 6.6424306340 5.2914044780 5.2914044780
+
+ 4.3860654351 4.3860654351
+
+ 4.3860654351
+
--- /dev/null
+4 6
+ 5.6053537261 6.2200032154 6.2159898496 6.4386968963
+ 6.2098101231 5.9596374798 5.4310222662 4.8790085257
+ 5.2638423822 5.4961954124 4.2928124162 4.3582250540
+ 4.2667755921 3.6813845958 3.5405429647 3.7026985912
+ 4.7177596156 3.3811992227 3.7936085458 4.5157490585
+ 6.2200032154 6.2159898496 5.2638423822 5.2638423822
+
+ 6.6296746680 6.6717260508 6.8168876945 6.8378960152
+ 6.3512039986 6.1425192423 5.5200206406 5.4858936024
+ 4.9382573262 4.2795129166 4.0534857018 4.2086428525
+ 3.5388086279 3.4209730327 2.7995744193 4.8164268553
+ 3.8790755071 3.5878262177 4.6457432555 6.6296746680
+ 6.6717260508 5.4858936024 5.4858936024
+
+ 6.6424306340 6.9715947250 6.9241225787 6.5291325353
+ 6.1394821460 5.4415971840 5.2914044780 4.7881880526
+ 4.1302408718 3.9405275117 3.9022731880 3.6078928145
+ 3.0809713206 3.0182595342 4.7935507272 3.8711928134
+ 3.6826168780 4.4872661030 6.6717260508 6.6424306340
+ 5.2914044780 5.2914044780
+
+ 7.3271707934 7.2970176793 7.0494948235 6.2565687662
+ 5.8133715113 5.8845382769 5.3463553241 4.7146928716
+ 4.4864793359 4.3092352779 3.6928361281 4.0100566871
+ 3.6584608489 4.5901676001 4.4681844010 4.6334334282
+ 4.9023945082 6.8168876945 6.9715947250 5.8845382769
+ 5.8845382769
+
+ 7.0106342493 6.8341607204 6.2047876082 5.5802122247
+ 5.6716968726 5.2011344932 4.1966522349 4.2383931955
+ 3.8954073453 3.5679086118 3.1136348127 2.7622561241
+ 4.5092053374 4.0238412606 3.9799853578 4.8087180377
+ 6.8378960152 6.9241225787 5.6716968726 5.6716968726
+
+ 6.5353606708 5.8844482033 5.1100976261 5.5515360893
+ 4.9587559238 4.2030790774 4.0351600709 3.6506289806
+ 3.4640416936 3.1611022087 2.7976192835 4.0422832763
+ 3.4415952740 3.8620694983 4.5934246009 6.3512039986
+ 6.5291325353 5.5515360893 5.5515360893
+
+ 5.2828811325 4.8298466315 4.7575327777 4.7249779042
+ 3.5656384285 3.5613366567 3.9273045564 3.7664779118
+ 3.5861699761 3.5523191627 4.6508672463 4.2424308260
+ 4.0987727783 4.5079669929 6.1425192423 6.1394821460
+ 4.7575327777 4.7575327777
+
+ 4.2222645754 4.1627149417 3.9786372270 3.1623319201
+ 2.9379696734 3.1365985188 2.8025851476 2.6576343487
+ 2.8733395837 3.8785891655 3.4443068252 3.4405132673
+ 4.0915305763 5.5200206406 5.4415971840 4.1627149417
+ 4.1627149417
+
+ 4.3860654351 3.6520014167 3.0264937881 2.5680579493
+ 2.5001857370 2.1393094782 1.4707077389 1.7373181916
+ 3.1045342752 2.0152329038 1.8196209136 3.6395852546
+ 5.4858936024 5.2914044780 4.3860654351 4.3860654351
+
+ 2.5016557935 2.7879606493 2.1424756816 1.0741544747
+ 1.3236438850 -.2712356678 .7613312181 2.8627348285
+ 1.5002127649 -.2392577636 3.5381128985 4.9382573262
+ 4.7881880526 3.6520014167 3.6520014167
+
+ 2.2903847447 2.2466771550 1.5053607856 1.4031797188
+ .9676209687 1.3075541842 2.8595118604 1.9833053538
+ -.0085773647 2.9557553666 4.2795129166 4.1302408718
+ 3.0264937881 3.0264937881
+
+ 1.2800038243 .7689806100 .9119606366 -.0744371793
+ .3669965427 2.4288563200 1.4975482174 -.2134101897
+ 2.9415021479 4.0534857018 3.9405275117 2.5680579493
+ 2.5680579493
+
+ -.6792428859 .4532383239 -.7590387660 -.3617034846
+ 1.6803275058 .6775209988 -.5354837468 2.6208591363
+ 4.2086428525 3.9022731880 2.5001857370 2.5001857370
+
+ .4136122500 -.1637152873 .2212564728 2.1287126020
+ .3896388275 -.2354089650 2.3255326077 3.5388086279
+ 3.6078928145 2.1393094782 2.1393094782
+
+ -3.5193531074 -2.0157307929 1.8355926160 2.7489174124
+ 2.0098391558 1.7975718667 3.4209730327 3.0809713206
+ 1.4707077389 1.4707077389
+
+ -1.3532688232 2.5504275249 2.8202873790 1.6253609671
+ 1.8625091247 2.7995744193 3.0182595342 1.7373181916
+ 1.7373181916
+
+ 3.7292778697 2.2944436481 -.0303565255 3.1115776177
+ 4.8164268553 4.7935507272 3.1045342752 3.1045342752
+
+ -.0827362961 -1.6043113182 2.4439837435 3.8790755071
+ 3.8711928134 2.0152329038 2.0152329038
+
+ -4.3897783280 2.3664634533 3.5878262177 3.6826168780
+ 1.8196209136 1.8196209136
+
+ 4.1927969260 4.6457432555 4.4872661030 3.6395852546
+ 3.6395852546
+
+ 6.6296746680 6.6717260508 5.4858936024 5.4858936024
+
+ 6.6424306340 5.2914044780 5.2914044780
+
+ 4.3860654351 4.3860654351
+
+ 4.3860654351
+
+ 2.6748060017 2.7338810145 2.9664647229 2.8819636737
+ 3.0210738150 2.8414286152 2.4773438660 2.4611943788
+ 2.4653201213 2.4925087371 2.5734767751 2.4564026744
+ 2.4847825281 2.4889289233 2.5089517645 2.5083338383
+ 2.4220622723 2.2714609770 2.4520703089 2.7026129788
+ 2.7338810145 2.9664647229 2.4653201213 2.4653201213
+
+ 4.9272154761 5.1054284230 4.2073516165 4.8513972837
+ 2.7848875293 3.5829861634 7.8660217576 7.4299209847
+ 1.9625939832 .7987769569 4.0580899681 1.8889021032
+ 3.1987197026 3.2673274538 2.6848131904 2.0043027404
+ 6.2446341910 8.1959452095 13.4748295858 2.6632376837
+ 5.1054284230 4.2073516165 1.9625939832 1.9625939832
+
+ .8699023011 1.0540660014 .9385909298 1.0263274101
+ 1.0835277045 1.0543183886 .7888686996 .8989305833
+ 1.0039962875 1.2427518128 .8932801724 .9173928990
+ 1.6157695657 1.4315860373 2.0498317879 1.4199615546
+ .9933677971 1.4319625600 27.4951763288 .7788025286
+ 1.0540660014 .9385909298 1.0039962875 1.0039962875
+
+ .0103697556 .0611385674 .0448303346 .0392831782
+ .0854166338 .0398896619 .0249496569 .0232410908
+ .0861379100 -.0754794185 -.0266146021 -.0163429099
+ .0572167103 -.0468608825 .0151048455 .0084963678
+ .0278930397 .0076922911 .1033536738 -.0098256036
+ .0611385674 .0448303346 .0861379100 .0861379100
+
+ 5.6053537261 6.2200032154 6.2159898496 6.4386968963
+ 6.2098101231 5.9596374798 5.4310222662 4.8790085257
+ 5.2638423822 5.4961954124 4.2928124162 4.3582250540
+ 4.2667755921 3.6813845958 3.5405429647 3.7026985912
+ 4.7177596156 3.3811992227 3.7936085458 4.5157490585
+ 6.2200032154 6.2159898496 5.2638423822 5.2638423822
+
+ 6.6296746680 6.6717260508 6.8168876945 6.8378960152
+ 6.3512039986 6.1425192423 5.5200206406 5.4858936024
+ 4.9382573262 4.2795129166 4.0534857018 4.2086428525
+ 3.5388086279 3.4209730327 2.7995744193 4.8164268553
+ 3.8790755071 3.5878262177 4.6457432555 6.6296746680
+ 6.6717260508 5.4858936024 5.4858936024
+
+ 6.6424306340 6.9715947250 6.9241225787 6.5291325353
+ 6.1394821460 5.4415971840 5.2914044780 4.7881880526
+ 4.1302408718 3.9405275117 3.9022731880 3.6078928145
+ 3.0809713206 3.0182595342 4.7935507272 3.8711928134
+ 3.6826168780 4.4872661030 6.6717260508 6.6424306340
+ 5.2914044780 5.2914044780
+
+ 7.3271707934 7.2970176793 7.0494948235 6.2565687662
+ 5.8133715113 5.8845382769 5.3463553241 4.7146928716
+ 4.4864793359 4.3092352779 3.6928361281 4.0100566871
+ 3.6584608489 4.5901676001 4.4681844010 4.6334334282
+ 4.9023945082 6.8168876945 6.9715947250 5.8845382769
+ 5.8845382769
+
+ 7.0106342493 6.8341607204 6.2047876082 5.5802122247
+ 5.6716968726 5.2011344932 4.1966522349 4.2383931955
+ 3.8954073453 3.5679086118 3.1136348127 2.7622561241
+ 4.5092053374 4.0238412606 3.9799853578 4.8087180377
+ 6.8378960152 6.9241225787 5.6716968726 5.6716968726
+
+ 6.5353606708 5.8844482033 5.1100976261 5.5515360893
+ 4.9587559238 4.2030790774 4.0351600709 3.6506289806
+ 3.4640416936 3.1611022087 2.7976192835 4.0422832763
+ 3.4415952740 3.8620694983 4.5934246009 6.3512039986
+ 6.5291325353 5.5515360893 5.5515360893
+
+ 5.2828811325 4.8298466315 4.7575327777 4.7249779042
+ 3.5656384285 3.5613366567 3.9273045564 3.7664779118
+ 3.5861699761 3.5523191627 4.6508672463 4.2424308260
+ 4.0987727783 4.5079669929 6.1425192423 6.1394821460
+ 4.7575327777 4.7575327777
+
+ 4.2222645754 4.1627149417 3.9786372270 3.1623319201
+ 2.9379696734 3.1365985188 2.8025851476 2.6576343487
+ 2.8733395837 3.8785891655 3.4443068252 3.4405132673
+ 4.0915305763 5.5200206406 5.4415971840 4.1627149417
+ 4.1627149417
+
+ 4.3860654351 3.6520014167 3.0264937881 2.5680579493
+ 2.5001857370 2.1393094782 1.4707077389 1.7373181916
+ 3.1045342752 2.0152329038 1.8196209136 3.6395852546
+ 5.4858936024 5.2914044780 4.3860654351 4.3860654351
+
+ 2.5016557935 2.7879606493 2.1424756816 1.0741544747
+ 1.3236438850 -.2712356678 .7613312181 2.8627348285
+ 1.5002127649 -.2392577636 3.5381128985 4.9382573262
+ 4.7881880526 3.6520014167 3.6520014167
+
+ 2.2903847447 2.2466771550 1.5053607856 1.4031797188
+ .9676209687 1.3075541842 2.8595118604 1.9833053538
+ -.0085773647 2.9557553666 4.2795129166 4.1302408718
+ 3.0264937881 3.0264937881
+
+ 1.2800038243 .7689806100 .9119606366 -.0744371793
+ .3669965427 2.4288563200 1.4975482174 -.2134101897
+ 2.9415021479 4.0534857018 3.9405275117 2.5680579493
+ 2.5680579493
+
+ -.6792428859 .4532383239 -.7590387660 -.3617034846
+ 1.6803275058 .6775209988 -.5354837468 2.6208591363
+ 4.2086428525 3.9022731880 2.5001857370 2.5001857370
+
+ .4136122500 -.1637152873 .2212564728 2.1287126020
+ .3896388275 -.2354089650 2.3255326077 3.5388086279
+ 3.6078928145 2.1393094782 2.1393094782
+
+ -3.5193531074 -2.0157307929 1.8355926160 2.7489174124
+ 2.0098391558 1.7975718667 3.4209730327 3.0809713206
+ 1.4707077389 1.4707077389
+
+ -1.3532688232 2.5504275249 2.8202873790 1.6253609671
+ 1.8625091247 2.7995744193 3.0182595342 1.7373181916
+ 1.7373181916
+
+ 3.7292778697 2.2944436481 -.0303565255 3.1115776177
+ 4.8164268553 4.7935507272 3.1045342752 3.1045342752
+
+ -.0827362961 -1.6043113182 2.4439837435 3.8790755071
+ 3.8711928134 2.0152329038 2.0152329038
+
+ -4.3897783280 2.3664634533 3.5878262177 3.6826168780
+ 1.8196209136 1.8196209136
+
+ 4.1927969260 4.6457432555 4.4872661030 3.6395852546
+ 3.6395852546
+
+ 6.6296746680 6.6717260508 5.4858936024 5.4858936024
+
+ 6.6424306340 5.2914044780 5.2914044780
+
+ 4.3860654351 4.3860654351
+
+ 4.3860654351
+
& faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,
& sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,
& eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,sigder,
- & dsci_inv,dscj_inv,gg
+ & dsci_inv,dscj_inv,gg,gg_lipi,gg_lipj
common /calc/ erij(3),rij,xj,yj,zj,dxi,dyi,dzi,dxj,dyj,dzj,
& chi1,chi2,chi12,chip1,chip2,chip12,alf1,alf2,alf12,om1,om2,om12,
& om1om2,chiom1,chiom2,chiom12,chipom1,chipom2,chipom12,eps1,
& faceps1,faceps1_inv,eps1_om12,facsig,sigsq,sigsq_om1,sigsq_om2,
& sigsq_om12,facp,facp_inv,facp1,eps2rt,eps2rt_om1,eps2rt_om2,
& eps2rt_om12,eps3rt,eom1,eom2,eom12,evdwij,eps2der,eps3der,sigder,
- & dsci_inv,dscj_inv,gg(3),i,j
+ & dsci_inv,dscj_inv,gg(3),gg_lipi(3),gg_lipj(3),i,j
& nsup,nstart_sup,nstart_seq,
& chain_length,iprzes,tabperm(maxperm,maxsym),nperm
common /from_zscore/ nz_start,nz_end,iz_sc
- double precision boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
- &buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
- common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad
- &buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+ double precision boxxsize,boxysize,boxzsize,enecut,sscut,
+ & sss,sssgrad,
+ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
+ common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad,
+ & buflipbot, bufliptop,bordlipbot,bordliptop,lipbufthick,lipthick
& gradcorr6_long,gcorr6_turn_long,gvdwx
integer nfl,icg
common /derivat/ dcdv(6,maxdim),dxdv(6,maxdim),dxds(6,maxres),
- & gradx(3,maxres,2),gradc(3,maxres,2),gvdwx(3,maxres),
- & gvdwc(3,maxres),gelc(3,maxres),gelc_long(3,maxres),
- & gvdwpp(3,maxres),gvdwc_scpp(3,maxres),gliptranc(3,maxres),
- & gliptranx(3,maxres),
- & gradx_scp(3,maxres),gvdwc_scp(3,maxres),ghpbx(3,maxres),
- & ghpbc(3,maxres),gloc(maxvar,2),gradcorr(3,maxres),
- & gradcorr_long(3,maxres),gradcorr5_long(3,maxres),
- & gradcorr6_long(3,maxres),gcorr6_turn_long(3,maxres),
- & gradxorr(3,maxres),gradcorr5(3,maxres),gradcorr6(3,maxres),
- & gloc_x(maxvar,2),gel_loc(3,maxres),gel_loc_long(3,maxres),
- & gcorr3_turn(3,maxres),
- &gcorr4_turn(3,maxres),gcorr6_turn(3,maxres),gradb(3,maxres),
- & gradbx(3,maxres),gel_loc_loc(maxvar),gel_loc_turn3(maxvar),
+ & gradx(3,-1:maxres,2),gradc(3,-1:maxres,2),gvdwx(3,-1:maxres),
+ & gvdwc(3,-1:maxres),gelc(3,-1:maxres),gelc_long(3,-1:maxres),
+ & gvdwpp(3,-1:maxres),gvdwc_scpp(3,-1:maxres),
+ & gliptranc(3,-1:maxres),
+ & gliptranx(3,-1:maxres),
+ & gradx_scp(3,-1:maxres),gvdwc_scp(3,-1:maxres),
+ & ghpbx(3,-1:maxres),
+ & ghpbc(3,-1:maxres),gloc(maxvar,2),gradcorr(3,-1:maxres),
+ & gradcorr_long(3,-1:maxres),gradcorr5_long(3,-1:maxres),
+ & gradcorr6_long(3,-1:maxres),gcorr6_turn_long(3,-1:maxres),
+ & gradxorr(3,-1:maxres),gradcorr5(3,-1:maxres),
+ & gradcorr6(3,-1:maxres),
+ & gloc_x(maxvar,2),gel_loc(3,-1:maxres),gel_loc_long(3,-1:maxres),
+ & gcorr3_turn(3,-1:maxres),
+ & gcorr4_turn(3,-1:maxres),gcorr6_turn(3,-1:maxres),
+ & gradb(3,-1:maxres),
+ & gradbx(3,-1:maxres),gel_loc_loc(maxvar),gel_loc_turn3(maxvar),
& gel_loc_turn4(maxvar),gel_loc_turn6(maxvar),
& gcorr_loc(maxvar),
- & g_corr5_loc(maxvar),g_corr6_loc(maxvar),gsccorc(3,maxres),
- & gsccorx(3,maxres),gsccor_loc(maxres),dtheta(3,2,maxres),
- & gscloc(3,maxres),gsclocx(3,maxres),
- & dphi(3,3,maxres),dalpha(3,3,maxres),domega(3,3,maxres),nfl,
+ & g_corr5_loc(maxvar),g_corr6_loc(maxvar),gsccorc(3,-1:maxres),
+ & gsccorx(3,-1:maxres),gsccor_loc(-1:maxres),
+ & dtheta(3,2,-1:maxres),
+ & gscloc(3,-1:maxres),gsclocx(3,-1:maxres),
+ & dphi(3,3,-1:maxres),dalpha(3,3,-1:maxres),domega(3,3,-1:maxres),
+ & nfl,
& icg
double precision derx,derx_turn
common /deriv_loc/ derx(3,5,2),derx_turn(3,5,2)
& iatel_s,iatel_e,iatscp_s,iatscp_e,iatel_s_vdw,iatel_e_vdw,
& ispp,iscp
common /interact/aa_aq(ntyp,ntyp),bb_aq(ntyp,ntyp),
- & aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp),aa(ntyp,ntyp),bb(ntyp,ntyp),
+ & aa_lip(ntyp,ntyp),bb_lip(ntyp,ntyp),
& augm(ntyp,ntyp),
& aad(ntyp,2),bad(ntyp,2),app(2,2),bpp(2,2),ael6(2,2),ael3(2,2),
& expon,expon2,nnt,nct,nint_gr(maxres),istart(maxres,maxint_gr),
& iatsc_s,iatsc_e,iatel_s,iatel_e,iatel_s_vdw,iatel_e_vdw,
& iatscp_s,iatscp_e,ispp,iscp
C 12/1/95 Array EPS included in the COMMON block.
- double precision eps,sigma,sigmaii,rs0,chi,chip,alp,sigma0,sigii,
+ double precision eps,epslip,sigma,sigmaii,rs0,chi,chip,alp,
+ & sigma0,sigii,
& rr0,r0,r0e,r0d,rpp,epp,elpp6,elpp3,eps_scp,rscp
- common /body/eps(ntyp,ntyp),sigma(0:ntyp1,0:ntyp1),
+ common /body/eps(ntyp,ntyp),epslip(ntyp,ntyp),
+ & sigma(0:ntyp1,0:ntyp1),
& sigmaii(ntyp,ntyp),
& rs0(ntyp,ntyp),chi(ntyp,ntyp),chip(ntyp),alp(ntyp),sigma0(ntyp),
& sigii(ntyp),rr0(ntyp),r0(ntyp,ntyp),r0e(ntyp,ntyp),r0d(ntyp,2),
include 'COMMON.HAIRPIN'
j=1
chain_beg=1
+ allareout=1
C do i=1,nres
C write(*,*) 'initial', i,j,c(j,i)
C enddo
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
evdw=evdw+(1.0d0-sss)*evdwij
C
rrij=1.0D0/rij
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
evdw=evdw+sss*evdwij
C
if (sss.lt.1.0d0) then
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
if (sss.gt.0.0d0) then
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*(1.0d0-sss)
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
C Calculate whole angle-dependent part of epsilon and contributions
C to its derivatives
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*(1.0d0-sss)
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=ssgradlipi*evdwij
+ gg_lipj(3)=ssgradlipj*evdwij
C Calculate angular part of the gradient.
call sc_grad_scale(1.0d0-sss)
endif
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=ssgradlipi*evdwij
+ gg_lipj(3)=ssgradlipj*evdwij
C Calculate angular part of the gradient.
call sc_grad_scale(sss)
endif
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+(evdwij+e_augm)*(1.0d0-sss)
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+(evdwij+e_augm)*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
enddo
c write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*scalfac
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*scalfac
c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
C Calculate the components of the gradient in DC and X
C
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
enddo
return
end
C
C Compute the side-chain and electrostatic interaction energy
C
+C print *,ipot
goto (101,102,103,104,105,106) ipot
C Lennard-Jones potential.
101 call elj(evdw)
goto 107
C Gay-Berne potential (shifted LJ, angular dependence).
104 call egb(evdw)
+C print *,"bylem w egb"
goto 107
C Gay-Berne-Vorobjev potential (shifted LJ, angular dependence).
105 call egbv(evdw)
#ifdef MPI
include 'mpif.h'
#endif
- double precision gradbufc(3,0:maxres),gradbufx(3,0:maxres),
- & glocbuf(4*maxres),gradbufc_sum(3,0:maxres),gloc_scbuf(3,0:maxres)
+ double precision gradbufc(3,-1:maxres),gradbufx(3,-1:maxres),
+ & glocbuf(4*maxres),gradbufc_sum(3,-1:maxres)
+ & ,gloc_scbuf(3,-1:maxres)
include 'COMMON.SETUP'
include 'COMMON.IOUNITS'
include 'COMMON.FFIELD'
call flush(iout)
#endif
#ifdef SPLITELE
- do i=1,nct
+ do i=0,nct
do j=1,3
gradbufc(j,i)=wsc*gvdwc(j,i)+
& wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
enddo
enddo
#else
- do i=1,nct
+ do i=0,nct
do j=1,3
gradbufc(j,i)=wsc*gvdwc(j,i)+
& wscp*(gvdwc_scp(j,i)+gvdwc_scpp(j,i))+
enddo
call flush(iout)
#endif
- do i=1,nres
+ do i=0,nres
do j=1,3
gradbufc_sum(j,i)=gradbufc(j,i)
enddo
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,nnt,-1
+ do i=nres-2,-1,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
enddo
call flush(iout)
#endif
- do i=1,nres
+ do i=-1,nres
do j=1,3
gradbufc_sum(j,i)=gradbufc(j,i)
gradbufc(j,i)=0.0d0
do j=1,3
gradbufc(j,nres-1)=gradbufc_sum(j,nres)
enddo
- do i=nres-2,nnt,-1
+ do i=nres-2,-1,-1
do j=1,3
gradbufc(j,i)=gradbufc(j,i+1)+gradbufc_sum(j,i+1)
enddo
do k=1,3
gradbufc(k,nres)=0.0d0
enddo
- do i=1,nct
+ do i=-1,nct
do j=1,3
#ifdef SPLITELE
gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+
eps0ij=eps(itypi,itypj)
fac=rrij**expon2
C have you changed here?
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
cd write (iout,'(2(a3,i3,2x),6(1pd12.4)/2(3(1pd12.4),5x)/)')
-cd & restyp(itypi),i,restyp(itypj),j,aa(itypi,itypj),
+cd & restyp(itypi),i,restyp(itypj),j,a(itypi,itypj),
cd & bb(itypi,itypj),1.0D0/dsqrt(rrij),evdwij,epsi,sigm,
cd & (c(k,i),k=1,3),(c(k,j),k=1,3)
evdw=evdw+evdwij
r_shift_inv=1.0D0/(rij+r0(itypi,itypj)-sigma(itypi,itypj))
fac=r_shift_inv**expon
C have you changed here?
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=e_augm+e1+e2
cd sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
cd epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
C to its derivatives
C have you changed here?
fac=(rrij*sigsq)**expon2
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),15(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
integer xshift,yshift,zshift
evdw=0.0D0
ccccc energy_dec=.false.
-c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
+C print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
lprn=.false.
c if (icall.eq.0) lprn=.false.
if (yi.lt.0) yi=yi+boxysize
zi=mod(zi,boxzsize)
if (zi.lt.0) zi=zi+boxzsize
+C define scaling factor for lipids
+
+C if (positi.le.0) positi=positi+boxzsize
+C print *,i
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
C xi=xi+xshift*boxxsize
C yi=yi+yshift*boxysize
C zi=zi+zshift*boxzsize
if (yj.lt.0) yj=yj+boxysize
zj=mod(zj,boxzsize)
if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((zj-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-zj)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)')
+ &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
+C if (ssgradlipj.gt.0.0d0) print *,"??WTF??"
+ print *,sslipi,sslipj,bordlipbot,zi,zj
dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
xj_safe=xj
yj_safe=yj
fac=rij_shift**expon
C here to start with
C if (c(i,3).gt.
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ faclip=fac
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
+C write(63,'(2i3,2e10.3,2f10.5)') i,j,aa,bb, evdwij,
+C &((sslipi+sslipj)/2.0d0+
+C &(2.0d0-sslipi-sslipj)/2.0d0)
c write (iout,*) "sigsq",sigsq," sig",sig," eps2rt",eps2rt,
c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij*sss
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,chi1,chi2,chip1,chip2,
fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
c fac=0.0d0
C Calculate the radial part of the gradient
+ gg_lipi(3)=eps1*(eps2rt*eps2rt)
+ &*(eps3rt*eps3rt)*sss/2.0d0*(faclip*faclip*
+ & (aa_lip(itypi,itypj)-aa_aq(itypi,itypj))
+ &+faclip*(bb_lip(itypi,itypj)-bb_aq(itypi,itypj)))
+ gg_lipj(3)=ssgradlipj*gg_lipi(3)
+ gg_lipi(3)=gg_lipi(3)*ssgradlipi
+C gg_lipi(3)=0.0d0
+C gg_lipj(3)=0.0d0
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+C define scaling factor for lipids
+
+C if (positi.le.0) positi=positi+boxzsize
+C print *,i
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+C xj=c(1,nres+j)-xi
+C yj=c(2,nres+j)-yi
+C zj=c(3,nres+j)-zi
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi-sslipj)/2.0d0
+C if (aa.ne.aa_aq(itypi,itypj)) write(63,'2e10.5')
+C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj))
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij+e_augm
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
write (iout,'(2(a3,i3,2x),17(0pf7.3))')
& restyp(itypi),i,restyp(itypj),j,
& epsi,sigm,sig,(augm(itypi,itypj)/epsi)**(1.0D0/12.0D0),
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac-2*expon*rrij*e_augm
+ fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
enddo
c write (iout,*) "gg",(gg(k),k=1,3)
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv*sss
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv*sss
c write (iout,*)(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
cgrad enddo
cgrad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(l)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(l)
enddo
return
end
include 'COMMON.CONTROL'
include 'COMMON.SPLITELE'
include 'COMMON.SBRIDGE'
+C this is done by Adasko
C print *,"wchodze"
C structure of box:
C water
C--bordlipbot--buffore ends
eliptran=0.0
do i=ilip_start,ilip_end
+C do i=1,1
if (itype(i).eq.ntyp1) cycle
- positi=(mod((c(3,i)+c(3,i+1)),boxzsize))
+ positi=(mod(((c(3,i)+c(3,i+1))/2.0d0),boxzsize))
if (positi.le.0) positi=positi+boxzsize
C print *,i
C first for peptide groups
sslip=sscalelip(fracinbuf)
ssgradlip=-sscagradlip(fracinbuf)/lipbufthick
eliptran=eliptran+sslip*pepliptran
- gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0
- gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0
+ gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
+ gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
C print *,"doing sccale for lower part"
gliptranc(3,i)=gliptranc(3,i)+ssgradlip*pepliptran/2.0d0
gliptranc(3,i-1)=gliptranc(3,i-1)+ssgradlip*pepliptran/2.0d0
C gliptranc(3,i-2)=gliptranc(3,i)+ssgradlip*pepliptran
- print *, "doing sscalefor top part"
+C print *, "doing sscalefor top part"
else
eliptran=eliptran+pepliptran
- print *,"I am in true lipid"
+C print *,"I am in true lipid"
endif
C else
C eliptran=elpitran+0.0 ! I am in water
C now multiply all by the peptide group transfer factor
C eliptran=eliptran*pepliptran
C now the same for side chains
+C do i=1,1
do i=ilip_start,ilip_end
if (itype(i).eq.ntyp1) cycle
positi=(mod(c(3,i+nres),boxzsize))
C print *,mod(c(3,i+nres),boxzsize),bordlipbot,bordliptop
c for each residue check if it is in lipid or lipid water border area
C respos=mod(c(3,i+nres),boxzsize)
+C print *,positi,bordlipbot,buflipbot
if ((positi.gt.bordlipbot)
& .and.(positi.lt.bordliptop)) then
C the energy transfer exist
eliptran=eliptran+sslip*liptranene(itype(i))
gliptranx(3,i)=gliptranx(3,i)
&+ssgradlip*liptranene(itype(i))/2.0d0
- gliptranc(3,i-1)=
- &+ssgradlip*liptranene(itype(i))
- print *,"doing sccale for lower part"
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))/2.0d0
+C print *,"doing sccale for lower part"
elseif (positi.gt.bufliptop) then
fracinbuf=1.0d0-
&((bordliptop-positi)/lipbufthick)
eliptran=eliptran+sslip*liptranene(itype(i))
gliptranx(3,i)=gliptranx(3,i)
&+ssgradlip*liptranene(itype(i))/2.0d0
- gliptranc(3,i-1)=
- &+ssgradlip*liptranene(itype(i))
- print *, "doing sscalefor top part",sslip,fracinbuf
+ gliptranc(3,i-1)= gliptranc(3,i-1)
+ &+ssgradlip*liptranene(itype(i))/2.0d0
+C print *, "doing sscalefor top part",sslip,fracinbuf
else
eliptran=eliptran+liptranene(itype(i))
- print *,"I am in true lipid"
+C print *,"I am in true lipid"
endif
endif ! if in lipid or buffor
C else
#ifdef DEBUG
write (iout,*) "gcart, gxcart, gloc before int_to_cart"
#endif
- do i=1,nct
+ do i=0,nct
do j=1,3
gcart(j,i)=gradc(j,i,icg)
gxcart(j,i)=gradx(j,i,icg)
C
C Initialize Cartesian-coordinate gradient
C
- do i=1,nres
+ do i=-1,nres
do j=1,3
gvdwx(j,i)=0.0D0
gradx_scp(j,i)=0.0D0
enddo
do i=1,ntyp
do j=1,ntyp
- aa(i,j)=0.0D0
- bb(i,j)=0.0D0
+ aa_aq(i,j)=0.0D0
+ bb_aq(i,j)=0.0D0
+ aa_lip(i,j)=0.0D0
+ bb_lip(i,j)=0.0D0
augm(i,j)=0.0D0
sigma(i,j)=0.0D0
r0(i,j)=0.0D0
include 'COMMON.MD'
include 'COMMON.IOUNITS'
include 'COMMON.SCCOR'
-c calculating dE/ddc1
+c calculating dE/ddc1
+C print *,"wchodze22",ialph(2,1)
if (nres.lt.3) go to 18
do j=1,3
gcart(j,1)=gcart(j,1)+gloc(1,icg)*dphi(j,1,4)
& gloc(ialph(2,1)+nside,icg)*domega(j,1,2)
endif
enddo
+C print *,"wchodze22",ialph(2,1)
c Calculating the remainder of dE/ddc2
do j=1,3
gcart(j,2)=gcart(j,2)+gloc(1,icg)*dphi(j,2,4)+
gcart(j,2)=gcart(j,2)+gloc(2,icg)*dphi(j,1,5)
endif
enddo
+C print *,"wchodze22",ialph(2,1)
c If there are only five residues
if(nres.eq.5) then
do j=1,3
endif
c If there are more than five residues
if(nres.gt.5) then
+C print *,"wchodze22",ialph(2,1)
do i=3,nres-3
+C print *,i,ialph(i,1)+nside
do j=1,3
gcart(j,i)=gcart(j,i)+gloc(i-2,icg)*dphi(j,3,i+1)
& +gloc(i-1,icg)*dphi(j,2,i+2)+
& gloc(i,icg)*dphi(j,1,i+3)+gloc(nres+i-4,icg)*dtheta(j,2,i+1)+
& gloc(nres+i-3,icg)*dtheta(j,1,i+2)
- if(itype(i).ne.10) then
+ if((itype(i).ne.10).and.(itype(i).ne.ntyp1)) then
gcart(j,i)=gcart(j,i)+gloc(ialph(i,1),icg)*dalpha(j,2,i)+
& gloc(ialph(i,1)+nside,icg)*domega(j,2,i)
endif
- if(itype(i+1).ne.10) then
+ if((itype(i+1).ne.10).and.(itype(i+1).ne.ntyp1)) then
gcart(j,i)=gcart(j,i)+gloc(ialph(i+1,1),icg)*dalpha(j,1,i+1)
& +gloc(ialph(i+1,1)+nside,icg)*domega(j,1,i+1)
endif
enddo
enddo
endif
+C print *,"wchodze22",ialph(2,1)
+
c Setting dE/ddnres-2
if(nres.gt.5) then
do j=1,3
c gloc(i,icg)=0.0D0
c write (iout,*) "poczotkoawy",i,gloc_sc(1,i,icg)
c enddo
+C print *,"tu dochodze??"
if (nres.lt.2) return
if ((nres.lt.3).and.(itype(1).eq.10)) return
if ((itype(1).ne.10).and.(itype(1).ne.ntyp1)) then
read (isidep,*,end=116,err=116)(sigii(i),i=1,ntyp)
read (isidep,*,end=116,err=116)(chip(i),i=1,ntyp)
read (isidep,*,end=116,err=116)(alp(i),i=1,ntyp)
+C now we start reading lipid
+ do i=1,ntyp
+ read (isidep,*,end=1161,err=1161)(epslip(i,j),j=i,ntyp)
+ print *,"WARNING!!"
+ do j=1,ntyp
+ epslip(i,j)=epslip(i,j)+0.5d0
+ enddo
+ enddo
C For the GB potential convert sigma'**2 into chi'
if (ipot.eq.4) then
do i=1,ntyp
epsij=eps(i,j)
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
- aa(i,j)=epsij*rrij*rrij
- bb(i,j)=-sigeps*epsij*rrij
- aa(j,i)=aa(i,j)
- bb(j,i)=bb(i,j)
+ aa_aq(i,j)=epsij*rrij*rrij
+ bb_aq(i,j)=-sigeps*epsij*rrij
+ aa_aq(j,i)=aa_aq(i,j)
+ bb_aq(j,i)=bb_aq(i,j)
+ epsijlip=epslip(i,j)
+ sigeps=dsign(1.0D0,epsijlip)
+ epsijlip=dabs(epsijlip)
+ aa_lip(i,j)=epsijlip*rrij*rrij
+ bb_lip(i,j)=-sigeps*epsijlip*rrij
+ aa_lip(j,i)=aa_lip(i,j)
+ bb_lip(j,i)=bb_lip(i,j)
if (ipot.gt.2) then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
endif
if (lprint) then
write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))')
- & restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
+ & restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),
& sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
endif
enddo
goto 999
116 write (iout,*) "Error reading electrostatic energy parameters."
goto 999
+ 1161 write (iout,*) "Error reading electrostatic energy parameters.Lip"
+ goto 999
117 write (iout,*) "Error reading side chain interaction parameters."
goto 999
118 write (iout,*) "Error reading SCp interaction parameters."
if ((lipbufthick*2.0d0).gt.lipthick)
&write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
endif
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+
+
if (me.eq.king .or. .not.out1file )
& write (iout,*) "DISTCHAINMAX",distchainmax
xi=c(1,nres+i)
yi=c(2,nres+i)
zi=c(3,nres+i)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
+
dxi=dc_norm(1,nres+i)
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
c alf1=0.0D0
c alf2=0.0D0
c alf12=0.0D0
- xj=c(1,nres+j)-xi
- yj=c(2,nres+j)-yi
- zj=c(3,nres+j)-zi
+C xj=c(1,nres+j)-xi
+C yj=c(2,nres+j)-yi
+C zj=c(3,nres+j)-zi
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij)
+ sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+ sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
+
C Calculate angle-dependent terms of energy and contributions to their
C derivatives.
call sc_angular
c---------------------------------------------------------------
rij_shift=1.0D0/rij_shift
fac=rij_shift**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
evdwij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=evdwij*eps3rt
eps3der=evdwij*eps2rt
evdwij=evdwij*eps2rt*eps3rt
evdw=evdw+evdwij
if (lprn) then
- sigm=dabs(aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
- epsi=bb(itypi,itypj)**2/aa(itypi,itypj)
+ sigm=dabs(aa/bb)**(1.0D0/6.0D0)
+ epsi=bb**2/aa
cd write (iout,'(2(a3,i3,2x),17(0pf7.3))')
cd & restyp(itypi),i,restyp(itypj),j,
cd & epsi,sigm,chi1,chi2,chip1,chip2,
fac=-expon*(e1+evdwij)*rij_shift
sigder=fac*sigder
fac=rij*fac
+ fac=fac+evdwij/sss*sssgrad/sigma(itypi,itypj)*rij
C Calculate the radial part of the gradient
gg(1)=xj*fac
gg(2)=yj*fac
gg(3)=zj*fac
+ gg_lipi(3)=ssgradlipi*evdwij
+ gg_lipj(3)=ssgradlipj*evdwij
C Calculate angular part of the gradient.
call sc_grad
ENDIF
dyi=dc_norm(2,nres+i)
dzi=dc_norm(3,nres+i)
dsci_inv=vbld_inv(i+nres)
-
+ xi=c(1,nres+i)
+ yi=c(2,nres+i)
+ zi=c(3,nres+i)
+ xi=mod(xi,boxxsize)
+ if (xi.lt.0) xi=xi+boxxsize
+ yi=mod(yi,boxysize)
+ if (yi.lt.0) yi=yi+boxysize
+ zi=mod(zi,boxzsize)
+ if (zi.lt.0) zi=zi+boxzsize
+C define scaling factor for lipids
+
+C if (positi.le.0) positi=positi+boxzsize
+C print *,i
+C first for peptide groups
+c for each residue check if it is in lipid or lipid water border area
+ if ((zi.gt.bordlipbot)
+ &.and.(zi.lt.bordliptop)) then
+C the energy transfer exist
+ if (zi.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipi=sscalelip(fracinbuf)
+ ssgradlipi=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipi=1.0d0
+ ssgradlipi=0.0
+ endif
+ else
+ sslipi=0.0d0
+ ssgradlipi=0.0
+ endif
itypj=itype(j)
- xj=c(1,nres+j)-c(1,nres+i)
- yj=c(2,nres+j)-c(2,nres+i)
- zj=c(3,nres+j)-c(3,nres+i)
+ xj=c(1,nres+j)
+ yj=c(2,nres+j)
+ zj=c(3,nres+j)
+ xj=mod(xj,boxxsize)
+ if (xj.lt.0) xj=xj+boxxsize
+ yj=mod(yj,boxysize)
+ if (yj.lt.0) yj=yj+boxysize
+ zj=mod(zj,boxzsize)
+ if (zj.lt.0) zj=zj+boxzsize
+ if ((zj.gt.bordlipbot)
+ &.and.(zj.lt.bordliptop)) then
+C the energy transfer exist
+ if (zj.lt.buflipbot) then
+C what fraction I am in
+ fracinbuf=1.0d0-
+ & ((positi-bordlipbot)/lipbufthick)
+C lipbufthick is thickenes of lipid buffore
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick
+ elseif (zi.gt.bufliptop) then
+ fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick)
+ sslipj=sscalelip(fracinbuf)
+ ssgradlipj=sscagradlip(fracinbuf)/lipbufthick
+ else
+ sslipj=1.0d0
+ ssgradlipj=0.0
+ endif
+ else
+ sslipj=0.0d0
+ ssgradlipj=0.0
+ endif
+ aa=aa_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +aa_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+ bb=bb_lip(itypi,itypj)*(sslipi+sslipj)/2.0d0
+ & +bb_aq(itypi,itypj)*(2.0d0-sslipi+sslipj)/2.0d0
+
+ dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ xj_safe=xj
+ yj_safe=yj
+ zj_safe=zj
+ subchap=0
+ do xshift=-1,1
+ do yshift=-1,1
+ do zshift=-1,1
+ xj=xj_safe+xshift*boxxsize
+ yj=yj_safe+yshift*boxysize
+ zj=zj_safe+zshift*boxzsize
+ dist_temp=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2
+ if(dist_temp.lt.dist_init) then
+ dist_init=dist_temp
+ xj_temp=xj
+ yj_temp=yj
+ zj_temp=zj
+ subchap=1
+ endif
+ enddo
+ enddo
+ enddo
+ if (subchap.eq.1) then
+ xj=xj_temp-xi
+ yj=yj_temp-yi
+ zj=zj_temp-zi
+ else
+ xj=xj_safe-xi
+ yj=yj_safe-yi
+ zj=zj_safe-zi
+ endif
+
+C xj=c(1,nres+j)-c(1,nres+i)
+C yj=c(2,nres+j)-c(2,nres+i)
+C zj=c(3,nres+j)-c(3,nres+i)
dxj=dc_norm(1,nres+j)
dyj=dc_norm(2,nres+j)
dzj=dc_norm(3,nres+j)
rrij=1.0D0/(xj*xj+yj*yj+zj*zj)
rij=dsqrt(rrij) ! sc_angular needs rij to really be the inverse
+ sss=sscale((1.0d0/rij)/sigma(itypi,itypj))
+ sssgrad=sscagrad((1.0d0/rij)/sigma(itypi,itypj))
c The following are set in sc_angular
c erij(1)=xj*rij
c erij(2)=yj*rij
ljXs=sig-sig0ij
ljA=eps1*eps2rt**2*eps3rt**2
- ljB=ljA*bb(itypi,itypj)
- ljA=ljA*aa(itypi,itypj)
- ljxm=ljXs+(-2.0D0*aa(itypi,itypj)/bb(itypi,itypj))**(1.0D0/6.0D0)
+ ljB=ljA*bb
+ ljA=ljA*aa
+ ljxm=ljXs+(-2.0D0*aa/bb)**(1.0D0/6.0D0)
ssXs=d0cm
deltat1=1.0d0-om1
c Stop and plot energy and derivative as a function of distance
if (checkstop) then
ssm=ssC-0.25D0*ssB*ssB/ssA
- ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
+ ljm=-0.25D0*ljB*bb/aa
if (ssm.lt.ljm .and.
& dabs(rij-0.5d0*(ssxm+ljxm)).lt.0.35d0*(ljxm-ssxm)) then
nicheck=1000
havebond=.false.
ljd=rij-ljXs
fac=(1.0D0/ljd)**expon
- e1=fac*fac*aa(itypi,itypj)
- e2=fac*bb(itypi,itypj)
+ e1=fac*fac*aa
+ e2=fac*bb
eij=eps1*eps2rt*eps3rt*(e1+e2)
eps2der=eij*eps3rt
eps3der=eij*eps2rt
- eij=eij*eps2rt*eps3rt
+ eij=eij*eps2rt*eps3rt*sss
sigder=-sig/sigsq
e1=e1*eps1*eps2rt**2*eps3rt**2
ed=-expon*(e1+eij)/ljd
sigder=ed*sigder
+ ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
eom1=eps2der*eps2rt_om1-2.0D0*alf1*eps3der+sigder*sigsq_om1
eom2=eps2der*eps2rt_om2+2.0D0*alf2*eps3der+sigder*sigsq_om2
eom12=eij*eps1_om12+eps2der*eps2rt_om12
havebond=.true.
ssd=rij-ssXs
eij=ssA*ssd*ssd+ssB*ssd+ssC
-
+ eij=eij*sss
ed=2*akcm*ssd+akct*deltat12
+ ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
pom1=akct*ssd
pom2=v1ss+2*v2ss*cosphi+3*v3ss*cosphi*cosphi
eom1=-2*akth*deltat1-pom1-om2*pom2
fac1=deltasq_inv*fac*(xm-rij)
fac2=deltasq_inv*fac*(rij-ssxm)
ed=delta_inv*(Ht*hd2-ssm*hd1)
+ eij=eij*sss
+ ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
eom1=fac1*d_ssxm(1)+fac2*d_xm(1)+h1*d_ssm(1)
eom2=fac1*d_ssxm(2)+fac2*d_xm(2)+h1*d_ssm(2)
eom12=fac1*d_ssxm(3)+fac2*d_xm(3)+h1*d_ssm(3)
else
havebond=.false.
- ljm=-0.25D0*ljB*bb(itypi,itypj)/aa(itypi,itypj)
- d_ljm(1)=-0.5D0*bb(itypi,itypj)/aa(itypi,itypj)*ljB
+ ljm=-0.25D0*ljB*bb/aa
+ d_ljm(1)=-0.5D0*bb/aa*ljB
d_ljm(2)=d_ljm(1)*(0.5D0*eps2rt_om2/eps2rt+alf2/eps3rt)
d_ljm(3)=d_ljm(1)*(0.5D0*eps1_om12+0.5D0*eps2rt_om12/eps2rt-
+ alf12/eps3rt)
fac1=deltasq_inv*fac*(ljxm-rij)
fac2=deltasq_inv*fac*(rij-xm)
ed=delta_inv*(ljm*hd2-Ht*hd1)
+ eij=eij*sss
+ ed=ed+eij/sss*sssgrad/sigma(itypi,itypj)*rij
eom1=fac1*d_xm(1)+fac2*d_ljxm(1)+h2*d_ljm(1)
eom2=fac1*d_xm(2)+fac2*d_ljxm(2)+h2*d_ljm(2)
eom12=fac1*d_xm(3)+fac2*d_ljxm(3)+h2*d_ljm(3)
checkstop=.false.
endif
c-------END TESTING CODE
+ gg_lipi(3)=ssgradlipi*eij
+ gg_lipj(3)=ssgradlipj*eij
do k=1,3
dcosom1(k)=(dc_norm(k,nres+i)-om1*erij(k))/rij
gg(k)=ed*erij(k)+eom1*dcosom1(k)+eom2*dcosom2(k)
enddo
do k=1,3
- gvdwx(k,i)=gvdwx(k,i)-gg(k)
+ gvdwx(k,i)=gvdwx(k,i)-gg(k)+gg_lipi(k)
& +(eom12*(dc_norm(k,nres+j)-om12*dc_norm(k,nres+i))
& +eom1*(erij(k)-om1*dc_norm(k,nres+i)))*dsci_inv
- gvdwx(k,j)=gvdwx(k,j)+gg(k)
+ gvdwx(k,j)=gvdwx(k,j)+gg(k)+gg_lipj(k)
& +(eom12*(dc_norm(k,nres+i)-om12*dc_norm(k,nres+j))
& +eom2*(erij(k)-om2*dc_norm(k,nres+j)))*dscj_inv
enddo
cgrad enddo
do l=1,3
- gvdwc(l,i)=gvdwc(l,i)-gg(l)
- gvdwc(l,j)=gvdwc(l,j)+gg(l)
+ gvdwc(l,i)=gvdwc(l,i)-gg(l)+gg_lipi(k)
+ gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k)
enddo
return
projects / unres.git / commitdiff