+! 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"
+
+! Ions by Aga
+
+ allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+ allocate(alphapolcat2(ntyp,ntyp))
+ allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+ allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+ allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
+ allocate(epsintabcat(ntyp,ntyp))
+ allocate(dtailcat(2,ntyp,ntyp))
+ allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
+ allocate(wqdipcat(2,ntyp,ntyp))
+ allocate(wstatecat(4,ntyp,ntyp))
+ allocate(dheadcat(2,2,ntyp,ntyp))
+ allocate(nstatecat(ntyp,ntyp))
+ allocate(debaykapcat(ntyp,ntyp))
+
+ if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
+ if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+
+ if (.not.allocated(ichargecat))&
+ allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
+ write(*,*) "before ions",oldion
+ ichargecat(:)=0
+
+! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
+ if (oldion.eq.0) then
+ if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
+ allocate(icharge(1:ntyp1))
+ read(iion,*) (icharge(i),i=1,ntyp)
+ else
+ read(iion,*) ijunk
+ endif
+ print *,ntyp_molec(5)
+ do i=-ntyp_molec(5),ntyp_molec(5)
+ read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
+ print *,msc(i,5),restok(i,5)
+ enddo
+ ip(5)=0.2
+ ! mp(5)=0.2
+ pstok(5)=3.0
+!DIR$ NOUNROLL
+ do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
+ do i=1,ntyp
+! do j=1,ntyp_molec(5)
+! write (*,*) "Im in ALAB", i, " ", j
+ read(iion,*) &
+ epscat(i,j),sigmacat(i,j), &
+! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
+ chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+
+ (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+! chiscat(i,j),chiscat(j,i), &
+ chis1cat(i,j),chis2cat(i,j), &
+
+ nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
+ dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+ dtailcat(1,i,j),dtailcat(2,i,j), &
+ epsheadcat(i,j),sig0headcat(i,j), &
+!wdipcat = w1 , w2
+! rborncat(i,j),rborncat(j,i),&
+ rborn1cat(i,j),rborn2cat(i,j),&
+ (wqdipcat(k,i,j),k=1,2), &
+ alphapolcat(i,j),alphapolcat2(j,i), &
+ (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
+
+ if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
+! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
+! if (i.eq.1) then
+! write (iout,*) 'i= ', i, ' j= ', j
+! write (iout,*) 'epsi0= ', epscat(1,j)
+! write (iout,*) 'sigma0= ', sigmacat(1,j)
+! write (iout,*) 'chi(i,j)= ', chicat(1,j)
+! write (iout,*) 'chi(j,i)= ', chicat(j,1)
+! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
+! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
+! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
+! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
+! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
+! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
+! write (iout,*) 'sig1= ', sigmap1cat(1,j)
+! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
+! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
+! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
+! write (iout,*) 'a1= ', rborncat(j,1)
+! write (iout,*) 'a2= ', rborncat(1,j)
+! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
+! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
+! write (iout,*) 'w1= ', wqdipcat(1,1,j)
+! write (iout,*) 'w2= ', wqdipcat(2,1,j)
+! endif
+
+!
+! If ((i.eq.1).and.(j.eq.27)) then
+! write (iout,*) 'SEP'
+! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+! endif
+
+ END DO
+ END DO
+ allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)
+ epsij=epscat(i,j)
+ rrij=sigmacat(i,j)
+ rrij=rrij**expon
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_aq_cat(i,j)=epsij*rrij*rrij
+ bb_aq_cat(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)-1
+ if (i.eq.10) then
+ write (iout,*) 'i= ', i, ' j= ', j
+ write (iout,*) 'epsi0= ', epscat(i,j)
+ write (iout,*) 'sigma0= ', sigmacat(i,j)
+ write (iout,*) 'chi1= ', chi1cat(i,j)
+ write (iout,*) 'chi1= ', chi2cat(i,j)
+ write (iout,*) 'chip1= ', chipp1cat(i,j)
+ write (iout,*) 'chip2= ', chipp2cat(i,j)
+ write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
+ write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
+ write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
+ write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
+ write (iout,*) 'sig1= ', sigmap1cat(i,j)
+ write (iout,*) 'sig2= ', sigmap2cat(i,j)
+ write (iout,*) 'chis1= ', chis1cat(i,j)
+ write (iout,*) 'chis1= ', chis2cat(i,j)
+ write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
+ write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
+ write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
+ write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
+ write (iout,*) 'a1= ', rborn1cat(i,j)
+ write (iout,*) 'a2= ', rborn2cat(i,j)
+ write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
+ write (iout,*) 'alphapol1= ', alphapolcat(i,j)
+ write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
+ write (iout,*) 'w1= ', wqdipcat(1,i,j)
+ write (iout,*) 'w2= ', wqdipcat(2,i,j)
+ write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
+ endif
+
+ If ((i.eq.1).and.(j.eq.27)) then
+ write (iout,*) 'SEP'
+ Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+ Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+ endif
+
+ enddo
+ enddo
+
+ endif
+! read number of Zn2+
+! here two denotes the Zn2+ and Cu2+
+ write(iout,*) "before TRANPARM"
+ allocate(aomicattr(0:3,2))
+ allocate(athetacattran(0:6,5,2))
+ allocate(agamacattran(3,5,2))
+ allocate(acatshiftdsc(5,2))
+ allocate(bcatshiftdsc(5,2))
+ allocate(demorsecat(5,2))
+ allocate(alphamorsecat(5,2))
+ allocate(x0catleft(5,2))
+ allocate(x0catright(5,2))
+ allocate(x0cattrans(5,2))
+ allocate(ntrantyp(2))
+ do i=1,1 ! currently only Zn2+
+
+ read(iiontran,*) ntrantyp(i)
+!ntrantyp=4
+!| ao0 ao1 ao2 ao3
+!ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+!CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
+!GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
+!HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
+ read(iiontran,*) (aomicattr(j,i),j=0,3)
+ do j=1,ntrantyp(i)
+ read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
+ (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
+ demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
+ x0cattrans(j,i)
+ enddo
+ enddo