+! MARTINI PARAMETER
+ allocate(ichargelipid(ntyp_molec(4)))
+ allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
+ allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
+ kjtokcal=0.2390057361
+ krad=57.295779513
+ !HERE THE MASS of MARTINI
+ write(*,*) "before MARTINI PARAM"
+ do i=1,ntyp_molec(4)
+ msc(i,4)=0.0d0
+ mp(4)=72.0d0
+ isc(i,4)=0.d0
+ enddo
+ ip(4)=0.0
+ !relative dielectric constant = 15 for implicit screening
+ k_coulomb_lip=332.0d0/15.0d0
+ !kbond = 1250 kJ/(mol*nm*2)
+ kbondlip=1250.0d0*kjtokcal/100.0d0
+ krad2=krad**2.0
+ lip_angle_force=0.0d0
+ if (DRY_MARTINI.gt.0) then
+ lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
+ lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
+ lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
+ lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
+ lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
+ lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
+ lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
+ lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
+ else
+ lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
+ lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
+ lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
+ lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
+ lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
+ lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
+ lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
+ lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
+ endif
+ lip_angle_angle=0.0d0
+ lip_angle_angle(3,12,12)=120.0/krad
+ lip_angle_angle(3,12,18)=180.0/krad
+ lip_angle_angle(3,18,16)=180.0/krad
+ lip_angle_angle(12,18,16)=180.0/krad
+ lip_angle_angle(18,16,18)=120.0/krad
+ lip_angle_angle(16,18,18)=180.0/krad
+ lip_angle_angle(12,18,18)=180.0/krad
+ lip_angle_angle(18,18,18)=180.0/krad
+ read(ilipbond,*) temp1
+ do i=1,18
+ read(ilipbond,*) temp1, lip_bond(i,1), &
+ lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
+ lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
+ lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
+ lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
+ lip_bond(i,18)
+ do j=1,18
+ lip_bond(i,j)=lip_bond(i,j)*10
+ enddo
+ enddo
+
+ read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
+ read(ilipnonbond,*) temp1
+ do i=1,18
+ read(ilipnonbond,*) temp1, lip_eps(i,1), &
+ lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
+ lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
+ lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
+ lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
+ lip_eps(i,18)
+! write(*,*) i, lip_eps(i,18)
+ do j=1,18
+ lip_eps(i,j)=lip_eps(i,j)*kjtokcal
+ enddo
+ enddo
+ read(ilipnonbond,*) temp1
+ do i=1,18
+ read(ilipnonbond,*) temp1,lip_sig(i,1), &
+ lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
+ lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
+ lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
+ lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
+ lip_sig(i,18)
+ do j=1,18
+ lip_sig(i,j)=lip_sig(i,j)*10.0
+ enddo
+ enddo
+ write(*,*) "after MARTINI PARAM"