From: Adam Sieradzan Date: Thu, 6 Aug 2015 07:01:16 +0000 (+0200) Subject: Merge branch 'AFM' into multichain X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=8df996bdc4bc56846b3799e9a2d47c24a9273125;hp=ad2687da4ea0c0ddde6c400a60c4817a689b3dff Merge branch 'AFM' into multichain --- diff --git a/CMakeLists.txt b/CMakeLists.txt index bb178aa..bd9465a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -195,7 +195,7 @@ if(UNRES_NA_MMCE) add_subdirectory(source/cluster/unres/src) add_subdirectory(source/xdrfpdb/src) add_subdirectory(source/xdrfpdb/src-M) - add_subdirectory(source/maxlik/src_CSA) +# add_subdirectory(source/maxlik/src_CSA) else() add_subdirectory(source/unres/src_MD) diff --git a/PARAM/Lip_tran_initial_ext.parm b/PARAM/Lip_tran_initial_ext.parm new file mode 100644 index 0000000..65acddb --- /dev/null +++ b/PARAM/Lip_tran_initial_ext.parm @@ -0,0 +1,25 @@ + 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 diff --git a/PARAM/scinter_GB_ext.parm b/PARAM/scinter_GB_ext.parm index 7f3e726..cae5b65 100644 --- a/PARAM/scinter_GB_ext.parm +++ b/PARAM/scinter_GB_ext.parm @@ -134,3 +134,112 @@ .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 + diff --git a/PARAM/scinter_GB_ext_lip.parm b/PARAM/scinter_GB_ext_lip.parm new file mode 100644 index 0000000..f4b894d --- /dev/null +++ b/PARAM/scinter_GB_ext_lip.parm @@ -0,0 +1,245 @@ +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 +C From here are lipid + 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 + diff --git a/source/unres/src_MD-M/COMMON.CALC b/source/unres/src_MD-M/COMMON.CALC index 67b4bb9..bf255c9 100644 --- a/source/unres/src_MD-M/COMMON.CALC +++ b/source/unres/src_MD-M/COMMON.CALC @@ -5,11 +5,11 @@ & 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 diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 777cc43..84ebf1b 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -1,7 +1,9 @@ integer nres,nsup,nstart_sup,nz_start,nz_end,iz_sc, - & nres0,nstart_seq,chain_length,iprzes,tabperm,nperm + & nres0,nstart_seq,chain_length,iprzes,tabperm,nperm,afmend, + & afmbeg double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r, - & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2 + & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2,velAFMconst, + & totTafm common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2), & dc_norm2(3,0:maxres2), @@ -14,5 +16,12 @@ & 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 - common /box/ boxxsize,boxysize,boxzsize,enecut,sscut,sss,sssgrad + 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 + common /afm/ forceAFMconst, distafminit,afmend,afmbeg, + & velAFMconst, + & totTafm + diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 40346c7..b8a775e 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -1,5 +1,5 @@ integer modecalc,iscode,indpdb,indback,indphi,iranconf,icheckgrad, - & inprint,i2ndstr,mucadyn,constr_dist,symetr + & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, @@ -8,6 +8,7 @@ & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, - & constr_dist,gnorm_check,gradout,split_ene,symetr + & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog, + & selfguide C... minim = .true. means DO minimization. C... energy_dec = .true. means print energy decomposition matrix diff --git a/source/unres/src_MD-M/COMMON.DERIV b/source/unres/src_MD-M/COMMON.DERIV index 7688aeb..8082b0a 100644 --- a/source/unres/src_MD-M/COMMON.DERIV +++ b/source/unres/src_MD-M/COMMON.DERIV @@ -1,27 +1,37 @@ - double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gelc_long, - & gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp, + double precision dcdv,dxdv,dxds,gradx,gradc,gvdwc,gelc,gelc_long + & gvdwpp,gel_loc,gel_loc_long,gvdwc_scpp,gliptranc,gliptranx, & gradx_scp,gvdwc_scp,ghpbx,ghpbc,gloc,gloc_x,dtheta,dphi,dalpha, & domega,gscloc,gsclocx,gradcorr,gradcorr_long,gradcorr5_long, & 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), - & gradx_scp(3,maxres),gvdwc_scp(3,maxres),ghpbx(3,maxres), - & ghpbc(3,maxres),gloc(0: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), - & 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,icg + & 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), + & gradafm(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,-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) double precision dXX_C1tab(3,maxres),dYY_C1tab(3,maxres), diff --git a/source/unres/src_MD-M/COMMON.FFIELD b/source/unres/src_MD-M/COMMON.FFIELD index d7d8cde..a63fe78 100644 --- a/source/unres/src_MD-M/COMMON.FFIELD +++ b/source/unres/src_MD-M/COMMON.FFIELD @@ -6,7 +6,7 @@ C----------------------------------------------------------------------- integer n_ene_comp,rescale_mode common /ffield/ wsc,wscp,welec,wbond,wstrain,wtor,wtor_d,wang, & wscloc,wcorr,wcorr4,wcorr5,wcorr6,wsccor,wel_loc,wturn3,wturn4, - & wturn6,wvdwpp,weights(n_ene),temp0, + & wturn6,wvdwpp,weights(n_ene),wliptran,temp0, & scal14,cutoff_corr,delt_corr,r0_corr,ipot,n_ene_comp, & rescale_mode common /potentials/ potname(5) diff --git a/source/unres/src_MD-M/COMMON.INTERACT b/source/unres/src_MD-M/COMMON.INTERACT index 37d3e88..448e829 100644 --- a/source/unres/src_MD-M/COMMON.INTERACT +++ b/source/unres/src_MD-M/COMMON.INTERACT @@ -1,11 +1,14 @@ - double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6 + double precision aa,bb,augm,aad,bad,app,bpp,ale6,ael3,ael6, + &aa_lip,bb_lip,aa_aq,bb_aq integer expon,expon2 integer nnt,nct,nint_gr,istart,iend,itype,itel,itypro, & ielstart,ielend,ielstart_vdw,ielend_vdw,nscp_gr,iscpstart, & iscpend,iatsc_s,iatsc_e, & iatel_s,iatel_e,iatscp_s,iatscp_e,iatel_s_vdw,iatel_e_vdw, & ispp,iscp - common /interact/aa(ntyp,ntyp),bb(ntyp,ntyp),augm(ntyp,ntyp), + common /interact/aa_aq(ntyp,ntyp),bb_aq(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), & iend(maxres,maxint_gr),itype(maxres),itel(maxres),itypro, @@ -15,9 +18,11 @@ & 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), @@ -31,3 +36,7 @@ c 12/5/03 modified 09/18/03 Bond stretching parameters. & vbldsc0(maxbondterm,ntyp),akp, & aksc(maxbondterm,ntyp),abond0(maxbondterm,ntyp), & distchainmax,nbondterm(ntyp) +C 01/29/15 Lipidic parameters + double precision pepliptran,liptranene + common /lipid/ pepliptran,liptranene(ntyp) + diff --git a/source/unres/src_MD-M/COMMON.IOUNITS b/source/unres/src_MD-M/COMMON.IOUNITS index 49b6db3..56e655e 100644 --- a/source/unres/src_MD-M/COMMON.IOUNITS +++ b/source/unres/src_MD-M/COMMON.IOUNITS @@ -39,10 +39,10 @@ C CSA I/O units & files C Parameter files character*256 bondname,thetname,rotname,torname,tordname, & fouriername,elename,sidename,scpname,sccorname,patname, - & thetname_pdb,rotname_pdb + & thetname_pdb,rotname_pdb,liptranname common /parfiles/ bondname,thetname,rotname,torname,tordname, & fouriername,elename,sidename,scpname,sccorname,patname, - & thetname_pdb,rotname_pdb + & thetname_pdb,rotname_pdb,liptranname character*3 pot C----------------------------------------------------------------------- C INP - main input file diff --git a/source/unres/src_MD-M/COMMON.LOCAL b/source/unres/src_MD-M/COMMON.LOCAL index 9f0627b..5d1ced7 100644 --- a/source/unres/src_MD-M/COMMON.LOCAL +++ b/source/unres/src_MD-M/COMMON.LOCAL @@ -40,7 +40,8 @@ C Virtual-bond lenghts & iphi_end,iphid_start,iphid_end,ibond_start,ibond_end, & ibondp_start,ibondp_end,ivec_start,ivec_end,iset_start,iset_end, & iturn3_start,iturn3_end,iturn4_start,iturn4_end,iint_start, - & iint_end,iphi1_start,iphi1_end,itau_start,itau_end, + & iint_end,iphi1_start,iphi1_end,itau_start,itau_end,ilip_start, + & ilip_end, & ibond_displ(0:max_fg_procs-1),ibond_count(0:max_fg_procs-1), & ithet_displ(0:max_fg_procs-1),ithet_count(0:max_fg_procs-1), & iphi_displ(0:max_fg_procs-1),iphi_count(0:max_fg_procs-1), @@ -56,6 +57,6 @@ C Virtual-bond lenghts & iint_end,iphi1_start,iphi1_end,iint_count,iint_displ,ivec_displ, & ivec_count,iset_displ,itau_start,itau_end, & iset_count,ibond_displ,ibond_count,ithet_displ,ithet_count, - & iphi_displ,iphi_count,iphi1_displ,iphi1_count + & iphi_displ,iphi_count,iphi1_displ,iphi1_count,ilip_start,ilip_end C Inverses of the actual virtual bond lengths common /invlen/ vbld_inv(maxres2) diff --git a/source/unres/src_MD-M/COMMON.SCCOR b/source/unres/src_MD-M/COMMON.SCCOR index 7952bd1..b3e6a6d 100644 --- a/source/unres/src_MD-M/COMMON.SCCOR +++ b/source/unres/src_MD-M/COMMON.SCCOR @@ -12,7 +12,7 @@ cc Parameters of the SCCOR term & nlor_sccor(-ntyp:ntyp,-ntyp:ntyp), & vlor1sccor(maxterm_sccor,20,20), & vlor2sccor(maxterm_sccor,20,20), - & vlor3sccor(maxterm_sccor,20,20),gloc_sc(3,0:maxres2,10), + & vlor3sccor(maxterm_sccor,20,20),gloc_sc(3,-1:maxres2,10), & dcostau(3,3,3,maxres2),dsintau(3,3,3,maxres2), & dtauangle(3,3,3,maxres2),dcosomicron(3,3,3,maxres2), & domicron(3,3,3,maxres2) diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index 84cf3fd..da975d3 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -6,17 +6,17 @@ ******************************************************************************** C Max. number of processors. integer maxprocs - parameter (maxprocs=2048) + parameter (maxprocs=1028) C Max. number of fine-grain processors integer max_fg_procs c parameter (max_fg_procs=maxprocs) - parameter (max_fg_procs=512) + parameter (max_fg_procs=256) C Max. number of coarse-grain processors integer max_cg_procs parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=800) + parameter (maxres=600) C Appr. max. number of interaction sites integer maxres2,maxres6,mmaxres2 parameter (maxres2=2*maxres,maxres6=6*maxres) @@ -95,7 +95,7 @@ C Max. number of conformations in the pool parameter (max_pool=10) C Number of energy components integer n_ene,n_ene2 - parameter (n_ene=21,n_ene2=2*n_ene) + parameter (n_ene=23,n_ene2=2*n_ene) C Number of threads in deformation integer max_thread,max_thread2 parameter (max_thread=4,max_thread2=2*max_thread) diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 6154487..dc58df8 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -196,7 +196,11 @@ c Variable time step algorithm. #endif endif if (ntwe.ne.0) then - if (mod(itime,ntwe).eq.0) call statout(itime) + if (mod(itime,ntwe).eq.0) then + call statout(itime) +C call enerprint(potEcomp) +C print *,itime,'AFM',Eafmforc,etot + endif #ifdef VOUT do j=1,3 v_work(j)=d_t(j,0) @@ -230,6 +234,8 @@ c Variable time step algorithm. #endif endif if (mod(itime,ntwx).eq.0) then + write(iout,*) 'time=',itime +C call check_ecartint call returnbox write (tytul,'("time",f8.2)') totT if(mdpdb) then @@ -514,6 +520,8 @@ c Second step of the velocity Verlet algorithm endif if (rattle) call rattle2 totT=totT+d_time + totTafm=totT +C print *,totTafm,"TU?" if (d_time.ne.d_time0) then d_time=d_time0 #ifndef LANG0 @@ -911,6 +919,7 @@ c Compute the complete potential energy potE=potEcomp(0)-potEcomp(20) c potE=energia_short(0)+energia_long(0) totT=totT+d_time + totTafm=totT c Calculate the kinetic and the total energy and the kinetic temperature call kinetic(EK) totE=EK+potE @@ -997,6 +1006,8 @@ c Applying velocity Verlet algorithm - step 1 to coordinates d_t(j,0)=d_t_old(j,0)+adt enddo do i=nnt,nct-1 +C SPYTAC ADAMA +C do i=0,nres do j=1,3 adt=d_a_old(j,i)*d_time adt2=0.5d0*adt @@ -1006,6 +1017,7 @@ c Applying velocity Verlet algorithm - step 1 to coordinates enddo enddo do i=nnt,nct +C do i=0,nres if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then inres=i+nres do j=1,3 @@ -1560,6 +1572,7 @@ c inquire(file=mremd_rst_name,exist=file_exist) endif call random_vel totT=0.0d0 + totTafm=totT endif else c Generate initial velocities @@ -1567,6 +1580,8 @@ c Generate initial velocities & write(iout,*) "Initial velocities randomly generated" call random_vel totT=0.0d0 +CtotTafm is the variable for AFM time which eclipsed during + totTafm=totT endif c rest2name = prefix(:ilen(prefix))//'.rst' if(me.eq.king.or..not.out1file)then @@ -1789,7 +1804,7 @@ c----------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.TIME1' - double precision xv,sigv,lowb,highb + double precision xv,sigv,lowb,highb,vec_afm(3) c Generate random velocities from Gaussian distribution of mean 0 and std of KT/m c First generate velocities in the eigenspace of the G matrix c write (iout,*) "Calling random_vel dimen dimen3",dimen,dimen3 @@ -1803,10 +1818,27 @@ c call flush(iout) lowb=-5*sigv highb=5*sigv d_t_work_new(ii)=anorm_distr(xv,sigv,lowb,highb) + c write (iout,*) "i",i," ii",ii," geigen",geigen(i), c & " d_t_work_new",d_t_work_new(ii) enddo enddo +C if (SELFGUIDE.gt.0) then +C distance=0.0 +C do j=1,3 +C vec_afm(j)=c(j,afmend)-c(j,afmbeg) +C distance=distance+vec_afm(j)**2 +C enddo +C distance=dsqrt(distance) +C do j=1,3 +C d_t_work_new(j+(afmbeg-1)*3)=-velAFMconst*vec_afm(j)/distance +C d_t_work_new(j+(afmend-1)*3)=velAFMconst*vec_afm(j)/distance +C write(iout,*) "myvel",d_t_work_new(j+(afmbeg-1)*3), +C & d_t_work_new(j+(afmend-1)*3) +C enddo + +C endif + c diagnostics c Ek1=0.0d0 c ii=0 diff --git a/source/unres/src_MD-M/MREMD.F b/source/unres/src_MD-M/MREMD.F index 5b0fc35..f29b6b5 100644 --- a/source/unres/src_MD-M/MREMD.F +++ b/source/unres/src_MD-M/MREMD.F @@ -1873,3 +1873,4 @@ c & (d_restart1(j,i+2*nres*il),j=1,3) if(me.eq.king) close(irest2) return end +c------------------------------------------ diff --git a/source/unres/src_MD-M/brown_step.F b/source/unres/src_MD-M/brown_step.F index 0be97f5..49652b8 100644 --- a/source/unres/src_MD-M/brown_step.F +++ b/source/unres/src_MD-M/brown_step.F @@ -381,6 +381,7 @@ c Calculate energy and forces potE=potEcomp(0)-potEcomp(20) call cartgrad totT=totT+d_time + totTafm=totT c Calculate the kinetic and total energy and the kinetic temperature call kinetic(EK) #ifdef MPI diff --git a/source/unres/src_MD-M/chainbuild.F b/source/unres/src_MD-M/chainbuild.F index 0721b72..157baa9 100644 --- a/source/unres/src_MD-M/chainbuild.F +++ b/source/unres/src_MD-M/chainbuild.F @@ -19,6 +19,10 @@ C Set lprn=.true. for debugging perbox=.false. fail=.false. print *, 'enter chainbuild' + call chainbuild_cart + return + end +#ifdef DEBUG if (perbox) then cost=dcos(theta(3)) sint=dsin(theta(3)) @@ -175,6 +179,7 @@ C endif return end +#endif c------------------------------------------------------------------------- subroutine orig_frame C diff --git a/source/unres/src_MD-M/checkder_p.F b/source/unres/src_MD-M/checkder_p.F index 32d2366..99f00bc 100644 --- a/source/unres/src_MD-M/checkder_p.F +++ b/source/unres/src_MD-M/checkder_p.F @@ -395,6 +395,7 @@ c write (iout,*) "etot11",etot11," etot12",etot12 c write(iout,'(2i5,2(a,f15.10))')i,j," etot",etot," etot1",etot1 dc(j,i)=ddc(j)-aincr call chainbuild_cart +C print *,c(j,i) c call int_from_cart1(.false.) if (.not.split_ene) then call etotal(energia1(0)) diff --git a/source/unres/src_MD-M/energy_p_new-sep_barrier.F b/source/unres/src_MD-M/energy_p_new-sep_barrier.F index 7d2d27f..c0b4c84 100644 --- a/source/unres/src_MD-M/energy_p_new-sep_barrier.F +++ b/source/unres/src_MD-M/energy_p_new-sep_barrier.F @@ -1,4 +1,33 @@ C----------------------------------------------------------------------- + double precision function sscalelip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscale=1.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscalelip=1.0d0+r*r*(2*r-3.0d0) +C else +C sscale=0d0 +C endif + return + end +C----------------------------------------------------------------------- + double precision function sscagradlip(r) + double precision r,gamm + include "COMMON.SPLITELE" +C if(r.lt.r_cut-rlamb) then +C sscagrad=0.0d0 +C else if(r.le.r_cut.and.r.ge.r_cut-rlamb) then +C gamm=(r-(r_cut-rlamb))/rlamb + sscagradlip=r*(6*r-6.0d0) +C else +C sscagrad=0.0d0 +C endif + return + end + +C----------------------------------------------------------------------- double precision function sscale(r) double precision r,gamm include "COMMON.SPLITELE" @@ -75,8 +104,8 @@ cd & 'iend=',iend(i,iint) 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 @@ -164,8 +193,8 @@ C Change 12/1/95 to calculate four-body interactions 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 @@ -248,8 +277,8 @@ 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) @@ -331,8 +360,8 @@ C 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) @@ -445,16 +474,16 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. 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, @@ -558,16 +587,16 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. 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, @@ -672,6 +701,33 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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 @@ -732,8 +788,8 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 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 @@ -742,8 +798,8 @@ c & " eps3rt",eps3rt," eps1",eps1," e1",e1," e2",e2 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, @@ -766,6 +822,8 @@ 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_scale(1.0d0-sss) endif @@ -854,6 +912,32 @@ c write (iout,*) "i",i," j", j," itype",itype(i),itype(j) 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 @@ -914,8 +998,8 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 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 @@ -924,8 +1008,8 @@ 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, @@ -948,6 +1032,8 @@ 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_scale(sss) endif @@ -1044,8 +1130,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! 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 @@ -1054,8 +1140,8 @@ c--------------------------------------------------------------- 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), @@ -1166,8 +1252,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! 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 @@ -1176,8 +1262,8 @@ c--------------------------------------------------------------- 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), @@ -1234,10 +1320,10 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 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)) @@ -1249,8 +1335,8 @@ C 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 diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 439fb81..3192819 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -99,6 +99,7 @@ c endif 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) @@ -112,6 +113,7 @@ C Berne-Pechukas potential (dilated LJ, angular dependence). 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) @@ -199,6 +201,7 @@ c print *,"Processor",myrank," computed UB" C C Calculate the SC local energy. C +C print *,"TU DOCHODZE?" call esc(escloc) c print *,"Processor",myrank," computed USC" C @@ -229,6 +232,7 @@ C else esccor=0.0d0 endif +C print *,"PRZED MULIt" c print *,"Processor",myrank," computed Usccorr" C C 12/1/95 Multi-body terms @@ -261,6 +265,19 @@ C after the equilibration time Uconst=0.0d0 Uconst_back=0.0d0 endif +C 01/27/2015 added by adasko +C the energy component below is energy transfer into lipid environment +C based on partition function +C print *,"przed lipidami" + if (wliptran.gt.0) then + call Eliptransfer(eliptran) + endif +C print *,"za lipidami" + if (AFMlog.gt.0) then + call AFMforce(Eafmforce) + else if (selfguide.gt.0) then + call AFMvel(Eafmforce) + endif #ifdef TIMING time_enecalc=time_enecalc+MPI_Wtime()-time00 #endif @@ -302,6 +319,8 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor + energia(22)=eliptran + energia(23)=Eafmforce c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" @@ -393,20 +412,23 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) + eliptran=energia(22) + Eafmforce=energia(23) #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran+Eafmforce #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+wliptran*eliptran + & +Eafmforce #endif energia(0)=etot c detecting NaNQ @@ -443,8 +465,9 @@ cMS$ATTRIBUTES C :: proc_proc #ifdef MPI include 'mpif.h' #endif - double precision gradbufc(3,maxres),gradbufx(3,maxres), - & glocbuf(4*maxres),gradbufc_sum(3,maxres),gloc_scbuf(3,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' @@ -497,7 +520,7 @@ c enddo 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))+ @@ -508,10 +531,13 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) + & +wliptran*gliptranc(j,i) + & +gradafm(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))+ @@ -523,6 +549,9 @@ c enddo & wcorr6*gradcorr6_long(j,i)+ & wturn6*gcorr6_turn_long(j,i)+ & wstrain*ghpbc(j,i) + & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + enddo enddo #endif @@ -536,7 +565,7 @@ c enddo 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 @@ -579,7 +608,7 @@ c 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 @@ -600,7 +629,7 @@ c 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 @@ -609,7 +638,7 @@ c 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 @@ -637,7 +666,7 @@ c 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)+ @@ -658,6 +687,8 @@ c enddo & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) + & +wliptran*gliptranc(j,i) + & +gradafm(j,i) #else gradc(j,i,icg)=gradbufc(j,i)+welec*gelc(j,i)+ & wel_loc*gel_loc(j,i)+ @@ -677,12 +708,16 @@ c enddo & wturn6*gcorr6_turn(j,i)+ & wsccor*gsccorc(j,i) & +wscloc*gscloc(j,i) + & +wliptran*gliptranc(j,i) + & +gradafm(j,i) + #endif gradx(j,i,icg)=wsc*gvdwx(j,i)+wscp*gradx_scp(j,i)+ & wbond*gradbx(j,i)+ & wstrain*ghpbx(j,i)+wcorr*gradxorr(j,i)+ & wsccor*gsccorx(j,i) & +wscloc*gsclocx(j,i) + & +wliptran*gliptranx(j,i) enddo enddo #ifdef DEBUG @@ -971,6 +1006,8 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) + eliptran=energia(22) + Eafmforce=energia(23) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, @@ -979,7 +1016,7 @@ C------------------------------------------------------------------------ & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, & edihcnstr,ebr*nss, - & Uconst,etot + & Uconst,eliptran,wliptran,Eafmforce,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1003,7 +1040,10 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETOT= ',1pE16.6,' (total)') + #else write (iout,10) evdw,wsc,evdw2,wscp,ees,welec, & estr,wbond,ebe,wang, @@ -1011,7 +1051,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, - & ebr*nss,Uconst,etot + & ebr*nss,Uconst,eliptran,wliptran,Eafmforc,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -1034,6 +1074,8 @@ C------------------------------------------------------------------------ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ + & 'ELT=',1pE16.6, ' WEIGHT=',1pD16.6,' (Lipid transfer energy)'/ + & 'EAFM= ',1pE16.6,' (atomic-force microscopy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif return @@ -1088,13 +1130,14 @@ C Change 12/1/95 to calculate four-body interactions c write (iout,*)'i=',i,' j=',j,' itypi=',itypi,' itypj=',itypj eps0ij=eps(itypi,itypj) fac=rrij**expon2 - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C have you changed here? + 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 @@ -1238,8 +1281,9 @@ C rij=1.0D0/r_inv_ij 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) +C have you changed here? + 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) @@ -1365,17 +1409,18 @@ C Calculate the angle-dependent terms of energy & contributions to derivatives. call sc_angular C Calculate whole angle-dependent part of epsilon and contributions 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, @@ -1425,7 +1470,7 @@ C 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. @@ -1473,6 +1518,35 @@ c endif 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 @@ -1557,6 +1631,36 @@ c 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- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.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 +C if (aa.ne.aa_aq(itypi,itypj)) write(63,'(2e10.5)') +C &(aa-aa_aq(itypi,itypj)),(bb-bb_aq(itypi,itypj)) +C if (ssgradlipj.gt.0.0d0) print *,"??WTF??" +C print *,sslipi,sslipj,bordlipbot,zi,zj dist_init=(xj-xi)**2+(yj-yi)**2+(zj-zi)**2 xj_safe=xj yj_safe=yj @@ -1625,18 +1729,24 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) c--------------------------------------------------------------- rij_shift=1.0D0/rij_shift fac=rij_shift**expon - e1=fac*fac*aa(itypi,itypj) - e2=fac*bb(itypi,itypj) +C here to start with +C if (c(i,3).gt. + 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, @@ -1658,12 +1768,20 @@ c & evdwij,fac,sigma(itypi,itypj),expon 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 C Calculate angular part of the gradient. call sc_grad - endif ! sss + endif ENDIF ! dyn_ss enddo ! j enddo ! iint @@ -1707,6 +1825,41 @@ c if (icall.eq.0) lprn=.true. 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- + & ((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 + dxi=dc_norm(1,nres+i) dyi=dc_norm(2,nres+i) dzi=dc_norm(3,nres+i) @@ -1743,9 +1896,74 @@ c chip12=0.0D0 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- + & ((zj-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslipj=sscalelip(fracinbuf) + ssgradlipj=-sscagradlip(fracinbuf)/lipbufthick + elseif (zj.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 +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) @@ -1766,8 +1984,8 @@ C I hate to put IF's in the loops, but here don't have another choice!!!! 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 @@ -1776,8 +1994,8 @@ c--------------------------------------------------------------- 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), @@ -1791,6 +2009,7 @@ C Calculate gradient components. 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 @@ -1901,10 +2120,10 @@ c write (iout,*) "eom1",eom1," eom2",eom2," eom12",eom12 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)) @@ -1921,8 +2140,8 @@ cgrad gvdwc(l,k)=gvdwc(l,k)+gg(l) 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 @@ -3063,6 +3282,8 @@ C Loop over i,i+2 and i,i+3 pairs of the peptide groups C C 14/01/2014 TURN3,TUNR4 does no go under periodic boundry condition do i=iturn3_start,iturn3_end + if (i.le.1) cycle +C write(iout,*) "tu jest i",i if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds & .or.((i+4).gt.nres) @@ -3097,6 +3318,7 @@ C end of changes by Ana num_cont_hb(i)=num_conti enddo do i=iturn4_start,iturn4_end + if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds & .or.((i+5).gt.nres) @@ -3164,6 +3386,7 @@ c c Loop over all pairs of interacting peptide groups except i,i+2 and i,i+3 c do i=iatel_s,iatel_e + if (i.le.1) cycle if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds & .or.((i+2).gt.nres) @@ -3220,7 +3443,8 @@ c endif c write (iout,*) 'i',i,' ielstart',ielstart(i),' ielend',ielend(i) num_conti=num_cont_hb(i) do j=ielstart(i),ielend(i) -c write (iout,*) i,j,itype(i),itype(j) +C write (iout,*) i,j + if (j.le.1) cycle if (itype(j).eq.ntyp1.or. itype(j+1).eq.ntyp1 C changes suggested by Ana to avoid out of bounds & .or.((j+2).gt.nres) @@ -4907,8 +5131,8 @@ c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj, c & dhpb(i),dhpb1(i),forcon(i) C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. - if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. - & iabs(itype(jjj)).eq.1) then +C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. +C & iabs(itype(jjj)).eq.1) then cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds if (.not.dyn_ss .and. i.le.nss) then @@ -4956,7 +5180,6 @@ cgrad enddo ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) enddo endif - endif enddo ehpb=0.5D0*ehpb return @@ -9919,4 +10142,199 @@ crc print *,((prod(i,j),i=1,2),j=1,2) return end +CCC---------------------------------------------- + subroutine Eliptransfer(eliptran) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + 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--bordliptop-- buffore starts +C--bufliptop--- here true lipid starts +C lipid +C--buflipbot--- lipid ends buffore starts +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))/2.0d0),boxzsize)) + 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 ((positi.gt.bordlipbot) + &.and.(positi.lt.bordliptop)) then +C the energy transfer exist + if (positi.lt.buflipbot) then +C what fraction I am in + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + 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" +C print *,i,sslip,fracinbuf,ssgradlip + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0-((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*pepliptran + 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 sscalefor top part" +C print *,i,sslip,fracinbuf,ssgradlip + else + eliptran=eliptran+pepliptran +C print *,"I am in true lipid" + endif +C else +C eliptran=elpitran+0.0 ! I am in water + endif + enddo +C print *, "nic nie bylo w lipidzie?" +C now multiply all by the peptide group transfer factor +C eliptran=eliptran*pepliptran +C now the same for side chains +CV do i=1,1 + do i=ilip_start,ilip_end + if (itype(i).eq.ntyp1) cycle + positi=(mod(c(3,i+nres),boxzsize)) + if (positi.le.0) positi=positi+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 + if (positi.lt.buflipbot) then + fracinbuf=1.0d0- + & ((positi-bordlipbot)/lipbufthick) +C lipbufthick is thickenes of lipid buffore + sslip=sscalelip(fracinbuf) + ssgradlip=-sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *,"doing sccale for lower part" + elseif (positi.gt.bufliptop) then + fracinbuf=1.0d0- + &((bordliptop-positi)/lipbufthick) + sslip=sscalelip(fracinbuf) + ssgradlip=sscagradlip(fracinbuf)/lipbufthick + eliptran=eliptran+sslip*liptranene(itype(i)) + gliptranx(3,i)=gliptranx(3,i) + &+ssgradlip*liptranene(itype(i)) + gliptranc(3,i-1)= gliptranc(3,i-1) + &+ssgradlip*liptranene(itype(i)) +C print *, "doing sscalefor top part",sslip,fracinbuf + else + eliptran=eliptran+liptranene(itype(i)) +C print *,"I am in true lipid" + endif + endif ! if in lipid or buffor +C else +C eliptran=elpitran+0.0 ! I am in water + enddo + return + end +C--------------------------------------------------------- +C AFM soubroutine for constant force + subroutine AFMforce(Eafmforce) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + real*8 diffafm(3) + dist=0.0d0 + Eafmforce=0.0d0 + do i=1,3 + diffafm(i)=c(i,afmend)-c(i,afmbeg) + dist=dist+diffafm(i)**2 + enddo + dist=dsqrt(dist) + Eafmforce=-forceAFMconst*(dist-distafminit) + do i=1,3 + gradafm(i,afmend-1)=-forceAFMconst*diffafm(i)/dist + gradafm(i,afmbeg-1)=forceAFMconst*diffafm(i)/dist + enddo +C print *,'AFM',Eafmforce + return + end +C--------------------------------------------------------- +C AFM subroutine with pseudoconstant velocity + subroutine AFMvel(Eafmforce) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.GEO' + include 'COMMON.VAR' + include 'COMMON.LOCAL' + include 'COMMON.CHAIN' + include 'COMMON.DERIV' + include 'COMMON.NAMES' + include 'COMMON.INTERACT' + include 'COMMON.IOUNITS' + include 'COMMON.CALC' + include 'COMMON.CONTROL' + include 'COMMON.SPLITELE' + include 'COMMON.SBRIDGE' + real*8 diffafm(3) +C Only for check grad COMMENT if not used for checkgrad +C totT=3.0d0 +C-------------------------------------------------------- +C print *,"wchodze" + dist=0.0d0 + Eafmforce=0.0d0 + do i=1,3 + diffafm(i)=c(i,afmend)-c(i,afmbeg) + dist=dist+diffafm(i)**2 + enddo + dist=dsqrt(dist) + Eafmforce=0.5d0*forceAFMconst + & *(distafminit+totTafm*velAFMconst-dist)**2 +C Eafmforce=-forceAFMconst*(dist-distafminit) + do i=1,3 + gradafm(i,afmend-1)=-forceAFMconst* + &(distafminit+totTafm*velAFMconst-dist) + &*diffafm(i)/dist + gradafm(i,afmbeg-1)=forceAFMconst* + &(distafminit+totTafm*velAFMconst-dist) + &*diffafm(i)/dist + enddo +C print *,'AFM',Eafmforce,totTafm*velAFMconst,dist + return + end diff --git a/source/unres/src_MD-M/geomout.F b/source/unres/src_MD-M/geomout.F index 266ad3d..0b711bd 100644 --- a/source/unres/src_MD-M/geomout.F +++ b/source/unres/src_MD-M/geomout.F @@ -445,6 +445,46 @@ c----------------------------------------------------------------- open(istat,file=statname,access="append") #endif #endif + if (AFMlog.gt.0) then + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)') + & itime,totT,EK,potE,totE, + & rms,frac,frac_nn,kinetic_T,t_bath,gyrate(), + & potEcomp(23),me + format1="a133" + else +C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,7f12.3,i5,$)') + & itime,totT,EK,potE,totE, + & kinetic_T,t_bath,gyrate(), + & potEcomp(23),me + format1="a114" + endif + else if (selfguide.gt.0) then + distance=0.0 + do j=1,3 + distance=distance+(c(j,afmend)-c(j,afmbeg))**2 + enddo + distance=dsqrt(distance) + if (refstr) then + call rms_nac_nnc(rms,frac,frac_nn,co,.false.) + write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, + & f9.2,i5,$)') + & itime,totT,EK,potE,totE, + & rms,frac,frac_nn,kinetic_T,t_bath,gyrate(), + & distance,potEcomp(23),me + format1="a133" +C print *,"CHUJOWO" + else +C print *,'A CHUJ',potEcomp(23) + write (line1,'(i10,f15.2,8f12.3,i5,$)') + & itime,totT,EK,potE,totE, + & kinetic_T,t_bath,gyrate(), + & distance,potEcomp(23),me + format1="a114" + endif + else if (refstr) then call rms_nac_nnc(rms,frac,frac_nn,co,.false.) write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') @@ -457,6 +497,7 @@ c----------------------------------------------------------------- & amax,kinetic_T,t_bath,gyrate(),me format1="a114" endif + endif if(usampl.and.totT.gt.eq_time) then write(line2,'(i5,2f9.4,300f7.4)') iset,uconst,uconst_back, & (qfrag(ii1),ii1=1,nfrag),(qpair(ii2),ii2=1,npair), diff --git a/source/unres/src_MD-M/gradient_p.F b/source/unres/src_MD-M/gradient_p.F index 2c670f2..acd4472 100644 --- a/source/unres/src_MD-M/gradient_p.F +++ b/source/unres/src_MD-M/gradient_p.F @@ -308,7 +308,7 @@ cd write(iout,*) 'calling int_to_cart' #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) @@ -349,7 +349,7 @@ C------------------------------------------------------------------------- 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 @@ -381,6 +381,9 @@ C gradx(j,i,icg)=0.0d0 gscloc(j,i)=0.0d0 gsclocx(j,i)=0.0d0 + gliptranc(j,i)=0.0d0 + gliptranx(j,i)=0.0d0 + gradafm(j,i)=0.0d0 do intertyp=1,3 gloc_sc(intertyp,i,icg)=0.0d0 enddo diff --git a/source/unres/src_MD-M/initialize_p.F b/source/unres/src_MD-M/initialize_p.F index 33b33a4..7ee3e42 100644 --- a/source/unres/src_MD-M/initialize_p.F +++ b/source/unres/src_MD-M/initialize_p.F @@ -119,6 +119,14 @@ C icsa_in=40 crc for ifc error 118 icsa_pdb=42 +C Lipidic input file for parameters range 60-79 + iliptranpar=60 +C input file for transfer sidechain and peptide group inside the +C lipidic environment if lipid is implicite + +C DNA input files for parameters range 80-99 +C Suger input files for parameters range 100-119 +C All-atom input files for parameters range 120-149 C C Set default weights of the energy terms. C @@ -146,8 +154,10 @@ c call memmon_print_usage() 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 @@ -246,20 +256,30 @@ C Initialize the bridge arrays jhpb(i)=0 enddo C Initialize correlation arrays - do i=-maxtor,maxtor + do i=1,maxres do k=1,2 b1(k,i)=0.0 b2(k,i)=0.0 b1tilde(k,i)=0.0 c b2tilde(k,i)=0.0 do j=1,2 +C CC(j,k,i)=0.0 +C Ctilde(j,k,i)=0.0 +C DD(j,k,i)=0.0 +C Dtilde(j,k,i)=0.0 + EE(j,k,i)=0.0 + enddo + enddo + enddo + do i=-maxtor,maxtor + do k=1,2 + do j=1,2 CC(j,k,i)=0.0 Ctilde(j,k,i)=0.0 DD(j,k,i)=0.0 Dtilde(j,k,i)=0.0 - EE(j,k,i)=0.0 enddo - enddo + enddo enddo C C Initialize timing. @@ -632,6 +652,8 @@ C Partition local interactions call int_bounds(nct-nnt,ibondp_start,ibondp_end) ibondp_start=ibondp_start+nnt ibondp_end=ibondp_end+nnt + call int_bounds(nres,ilip_start,ilip_end) + ilip_start=ilip_start call int_bounds1(nres-1,ivec_start,ivec_end) c print *,"Processor",myrank,fg_rank,fg_rank1, c & " ivec_start",ivec_start," ivec_end",ivec_end @@ -1157,6 +1179,8 @@ C ibondp_end=nct-1 iset_end=nres+1 iint_start=2 iint_end=nres-1 + ilip_start=1 + ilip_end=nres #endif return end diff --git a/source/unres/src_MD-M/int_to_cart.f b/source/unres/src_MD-M/int_to_cart.f index f413622..d3a8a92 100644 --- a/source/unres/src_MD-M/int_to_cart.f +++ b/source/unres/src_MD-M/int_to_cart.f @@ -14,7 +14,9 @@ c------------------------------------------------------------- include 'COMMON.MD' include 'COMMON.IOUNITS' include 'COMMON.SCCOR' -c calculating dE/ddc1 + include 'COMMON.CONTROL' +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) @@ -24,6 +26,7 @@ c calculating dE/ddc1 & 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)+ @@ -40,6 +43,7 @@ c Calculating the remainder of dE/ddc2 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 @@ -58,23 +62,27 @@ c If there are only five residues 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 @@ -123,6 +131,7 @@ c do i=1,nres 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 @@ -267,6 +276,13 @@ c Settind dE/ddnres enddo endif c The side-chain vector derivatives +C if (SELFGUIDE.gt.0) then +C do j=1,3 +C gcart(j,afmbeg)=gcart(j,afmbeg)+gcart(j,afmend) +C gcart(j,afmbeg)=0.0d0 +C gcart(j,afmend)=0.0d0 +C enddo +C endif return end diff --git a/source/unres/src_MD-M/parmread.F b/source/unres/src_MD-M/parmread.F index 43ceeac..0620acb 100644 --- a/source/unres/src_MD-M/parmread.F +++ b/source/unres/src_MD-M/parmread.F @@ -97,6 +97,12 @@ c enddo enddo endif +C reading lipid parameters + read(iliptranpar,*) pepliptran + do i=1,ntyp + read(iliptranpar,*) liptranene(i) + enddo + close(iliptranpar) #ifdef CRYST_THETA C C Read the parameters of the probability distribution/energy expression @@ -1000,6 +1006,7 @@ c Dtilde(2,2,i)=0.0d0 EEold(2,1,-i)=-b(12,i)+b(13,i) EEold(1,2,-i)=-b(12,i)-b(13,i) write(iout,*) "TU DOCHODZE" + print *,"JESTEM" c ee(1,1,i)=1.0d0 c ee(2,2,i)=1.0d0 c ee(2,1,i)=0.0d0 @@ -1107,6 +1114,14 @@ C---------------------- GB or BP potential ----------------------------- 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.05d0 + enddo + enddo C For the GB potential convert sigma'**2 into chi' if (ipot.eq.4) then do i=1,ntyp @@ -1145,6 +1160,7 @@ C Calculate the "working" parameters of SC interactions. do i=2,ntyp do j=1,i-1 eps(i,j)=eps(j,i) + epslip(i,j)=epslip(j,i) enddo enddo do i=1,ntyp @@ -1173,10 +1189,17 @@ C Calculate the "working" parameters of SC interactions. 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 @@ -1209,7 +1232,7 @@ c augm(i,j)=0.5D0**(2*expon)*aa(i,j) 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 @@ -1313,6 +1336,8 @@ c v3ss=0.0d0 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." diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index 5da38bd..2026425 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -80,7 +80,7 @@ C Calculate the CM of the preceding residue. endif C Start new residue. c write (iout,'(a80)') card - read (card(24:26),*) ires + read (card(23:26),*) ires read (card(18:20),'(a3)') res if (ibeg.eq.1) then ishift=ires-1 @@ -369,7 +369,6 @@ cc enddiag hfrag(i,j)=hfrag(i,j)-ishift enddo enddo - return end c--------------------------------------------------------------------------- diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 68d63c4..6ab2a4e 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -137,6 +137,9 @@ C Set up the time limit (caution! The time must be input in minutes!) refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0) indpdb=index(controlcard,'PDBSTART') extconf=(index(controlcard,'EXTCONF').gt.0) + AFMlog=(index(controlcard,'AFM')) + selfguide=(index(controlcard,'SELFGUIDE')) + print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'IPRINT',iprint,0) call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) @@ -224,6 +227,25 @@ C Reading the dimensions of box in x,y,z coordinates c Cutoff range for interactions call reada(controlcard,"R_CUT",r_cut,15.0d0) call reada(controlcard,"LAMBDA",rlamb,0.3d0) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick +C endif + if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0)) + & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + 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 @@ -586,6 +608,7 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) + call reada(weightcard,'WLT',wliptran,0.0D0) if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WCORRH',wcorr,1.0D0) @@ -935,6 +958,7 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb + if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp call hpb_partition if(me.eq.king.or..not.out1file) & write (iout,*) 'Contact order:',co @@ -1972,6 +1996,8 @@ C Get parameter filenames and open the parameter files. open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) open (irotam_pdb,file=rotname_pdb,status='old',action='read') @@ -2176,6 +2202,7 @@ c------------------------------------------------------------------------------- include 'COMMON.MD' open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath + totTafm=totT do i=1,2*nres read(irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo @@ -2231,6 +2258,36 @@ c------------------------------------------------------------------------------- enddo return end +C--------------------------------------------------------------------------- + subroutine read_afminp + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + character*320 afmcard + print *, "wchodze" + call card_concat(afmcard) + call readi(afmcard,"BEG",afmbeg,0) + call readi(afmcard,"END",afmend,0) + call reada(afmcard,"FORCE",forceAFMconst,0.0d0) + call reada(afmcard,"VEL",velAFMconst,0.0d0) + print *,'FORCE=' ,forceAFMconst +CCCC NOW PROPERTIES FOR AFM + distafminit=0.0d0 + do i=1,3 + distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit + enddo + distafminit=dsqrt(distafminit) + print *,'initdist',distafminit + return + end + c------------------------------------------------------------------------------- subroutine read_dist_constr implicit real*8 (a-h,o-z) diff --git a/source/unres/src_MD-M/sc_move.F b/source/unres/src_MD-M/sc_move.F index 274767b..c552ee0 100644 --- a/source/unres/src_MD-M/sc_move.F +++ b/source/unres/src_MD-M/sc_move.F @@ -724,6 +724,35 @@ c if (icall.eq.0) lprn=.true. 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) @@ -757,14 +786,85 @@ c chip12=0.0D0 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 @@ -783,16 +883,16 @@ cd & rij_shift,1.0D0/rij,sig,sig0ij,sigsq,1-dsqrt(sigsq) 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, @@ -809,10 +909,13 @@ C Calculate gradient components. 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 diff --git a/source/unres/src_MD-M/ssMD.F b/source/unres/src_MD-M/ssMD.F index eab3c70..3c89ef5 100644 --- a/source/unres/src_MD-M/ssMD.F +++ b/source/unres/src_MD-M/ssMD.F @@ -150,11 +150,119 @@ c-------END TESTING CODE 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) @@ -172,6 +280,8 @@ c-------END TESTING CODE 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 @@ -187,9 +297,9 @@ c om12=dxi*dxj+dyi*dyj+dzi*dzj 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 @@ -223,7 +333,7 @@ c-------TESTING CODE 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 @@ -248,17 +358,18 @@ c-------END TESTING CODE 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 @@ -267,8 +378,9 @@ c-------END TESTING CODE 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 @@ -309,13 +421,15 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE 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) @@ -331,6 +445,8 @@ c-------FIRST METHOD, DISCONTINUOUS SECOND DERIVATIVE 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) @@ -433,6 +549,8 @@ c-------TESTING CODE 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 @@ -442,10 +560,10 @@ c-------END TESTING CODE 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 @@ -456,8 +574,8 @@ 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(k) + gvdwc(l,j)=gvdwc(l,j)+gg(l)+gg_lipj(k) enddo return