8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
748 sigmasctube,krad2,ract
749 integer :: ichir1,ichir2,ijunk,irdiff
751 character*80 temp1,mychar
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
817 allocate(isc(-ntyp-1:ntyp+1,5)) !(ntyp+1)
818 allocate(restok(-ntyp-1:ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
873 if (.not.allocated(ichargecat)) &
874 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
876 if (oldion.eq.1) then
878 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
879 print *,msc(i,5),restok(i,5)
883 read (iion,*) ncatprotparm
884 allocate(catprm(ncatprotparm,4))
886 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
890 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
894 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
897 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
898 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
902 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
903 restyp(i,5),"-",restyp(j,2)
906 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
907 !----------------------------------------------------
908 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
909 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
910 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
911 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
912 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
913 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
923 allocate(liptranene(ntyp))
924 !C reading lipid parameters
925 write (iout,*) "iliptranpar",iliptranpar
927 read(iliptranpar,*) pepliptran
930 read(iliptranpar,*) liptranene(i)
931 print *,liptranene(i)
934 ! water parmaters entalphy
935 allocate(awaterenta(0:400))
936 allocate(bwaterenta(0:400))
937 allocate(cwaterenta(0:400))
938 allocate(dwaterenta(0:400))
939 allocate(awaterentro(0:400))
940 allocate(bwaterentro(0:400))
941 allocate(cwaterentro(0:400))
942 allocate(dwaterentro(0:400))
944 read(iwaterwater,*) mychar
945 read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
946 cwaterenta(0),dwaterenta(0)
948 read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
949 cwaterenta(i),dwaterenta(i)
950 irdiff=int((ract-2.06d0)*50.0d0)+1
951 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
953 ! water parmaters entrophy
954 read(iwaterwater,*) mychar
955 read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
956 cwaterentro(0),dwaterentro(0)
958 read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
959 cwaterentro(i),dwaterentro(i)
960 irdiff=int((ract-2.06d0)*50.0d0)+1
961 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
967 ! Read the parameters of the probability distribution/energy expression
968 ! of the virtual-bond valence angles theta
971 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
972 (bthet(j,i,1,1),j=1,2)
973 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
974 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
975 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
979 athet(1,i,1,-1)=athet(1,i,1,1)
980 athet(2,i,1,-1)=athet(2,i,1,1)
981 bthet(1,i,1,-1)=-bthet(1,i,1,1)
982 bthet(2,i,1,-1)=-bthet(2,i,1,1)
983 athet(1,i,-1,1)=-athet(1,i,1,1)
984 athet(2,i,-1,1)=-athet(2,i,1,1)
985 bthet(1,i,-1,1)=bthet(1,i,1,1)
986 bthet(2,i,-1,1)=bthet(2,i,1,1)
990 athet(1,i,-1,-1)=athet(1,-i,1,1)
991 athet(2,i,-1,-1)=-athet(2,-i,1,1)
992 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
993 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
994 athet(1,i,-1,1)=athet(1,-i,1,1)
995 athet(2,i,-1,1)=-athet(2,-i,1,1)
996 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
997 bthet(2,i,-1,1)=bthet(2,-i,1,1)
998 athet(1,i,1,-1)=-athet(1,-i,1,1)
999 athet(2,i,1,-1)=athet(2,-i,1,1)
1000 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1001 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1002 theta0(i)=theta0(-i)
1006 polthet(j,i)=polthet(j,-i)
1009 gthet(j,i)=gthet(j,-i)
1015 if (.not.LaTeX) then
1016 write (iout,'(a)') &
1017 'Parameters of the virtual-bond valence angles:'
1018 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
1019 ' ATHETA0 ',' A1 ',' A2 ',&
1022 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1023 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
1025 write (iout,'(/a/9x,5a/79(1h-))') &
1026 'Parameters of the expression for sigma(theta_c):',&
1027 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
1028 ' ALPH3 ',' SIGMA0C '
1030 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1031 (polthet(j,i),j=0,3),sigc0(i)
1033 write (iout,'(/a/9x,5a/79(1h-))') &
1034 'Parameters of the second gaussian:',&
1035 ' THETA0 ',' SIGMA0 ',' G1 ',&
1038 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1039 sig0(i),(gthet(j,i),j=1,3)
1042 write (iout,'(a)') &
1043 'Parameters of the virtual-bond valence angles:'
1044 write (iout,'(/a/9x,5a/79(1h-))') &
1045 'Coefficients of expansion',&
1046 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1047 ' b1*10^1 ',' b2*10^1 '
1049 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1050 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1051 (10*bthet(j,i,1,1),j=1,2)
1053 write (iout,'(/a/9x,5a/79(1h-))') &
1054 'Parameters of the expression for sigma(theta_c):',&
1055 ' alpha0 ',' alph1 ',' alph2 ',&
1056 ' alhp3 ',' sigma0c '
1058 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1059 (polthet(j,i),j=0,3),sigc0(i)
1061 write (iout,'(/a/9x,5a/79(1h-))') &
1062 'Parameters of the second gaussian:',&
1063 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1066 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1067 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1073 ! Read the parameters of Utheta determined from ab initio surfaces
1074 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1076 IF (tor_mode.eq.0) THEN
1077 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1078 ntheterm3,nsingle,ndouble
1079 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1081 !----------------------------------------------------
1082 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1083 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1084 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1085 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1086 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1087 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1088 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1089 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1090 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1091 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1092 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1093 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1094 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1095 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1096 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1097 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1098 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1099 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1100 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1101 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1102 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1103 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1104 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1105 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1107 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1109 ithetyp(i)=-ithetyp(-i)
1112 aa0thet(:,:,:,:)=0.0d0
1113 aathet(:,:,:,:,:)=0.0d0
1114 bbthet(:,:,:,:,:,:)=0.0d0
1115 ccthet(:,:,:,:,:,:)=0.0d0
1116 ddthet(:,:,:,:,:,:)=0.0d0
1117 eethet(:,:,:,:,:,:)=0.0d0
1118 ffthet(:,:,:,:,:,:,:)=0.0d0
1119 ggthet(:,:,:,:,:,:,:)=0.0d0
1121 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1123 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1124 ! VAR:1=non-glicyne non-proline 2=proline
1125 ! VAR:negative values for D-aminoacid
1127 do j=-nthetyp,nthetyp
1128 do k=-nthetyp,nthetyp
1129 read (ithep,'(6a)',end=111,err=111) res1
1130 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1131 ! VAR: aa0thet is variable describing the average value of Foureir
1132 ! VAR: expansion series
1133 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1134 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1135 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1136 read (ithep,*,end=111,err=111) &
1137 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1138 read (ithep,*,end=111,err=111) &
1139 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1140 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1141 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1142 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1144 read (ithep,*,end=111,err=111) &
1145 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1146 ffthet(lll,llll,ll,i,j,k,iblock),&
1147 ggthet(llll,lll,ll,i,j,k,iblock),&
1148 ggthet(lll,llll,ll,i,j,k,iblock),&
1149 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1154 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1155 ! coefficients of theta-and-gamma-dependent terms are zero.
1156 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1157 ! RECOMENTDED AFTER VERSION 3.3)
1161 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1162 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1164 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1165 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1168 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1170 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1173 ! AND COMMENT THE LOOPS BELOW
1177 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1178 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1180 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1181 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1184 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1186 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1191 ! Substitution for D aminoacids from symmetry.
1194 do j=-nthetyp,nthetyp
1195 do k=-nthetyp,nthetyp
1196 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1198 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1202 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1203 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1204 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1205 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1211 ffthet(llll,lll,ll,i,j,k,iblock)= &
1212 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1213 ffthet(lll,llll,ll,i,j,k,iblock)= &
1214 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1215 ggthet(llll,lll,ll,i,j,k,iblock)= &
1216 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1217 ggthet(lll,llll,ll,i,j,k,iblock)= &
1218 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1227 ! Control printout of the coefficients of virtual-bond-angle potentials
1230 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1235 write (iout,'(//4a)') &
1236 'Type ',onelett(i),onelett(j),onelett(k)
1237 write (iout,'(//a,10x,a)') " l","a[l]"
1238 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1239 write (iout,'(i2,1pe15.5)') &
1240 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1242 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1243 "b",l,"c",l,"d",l,"e",l
1245 write (iout,'(i2,4(1pe15.5))') m,&
1246 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1247 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1251 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1252 "f+",l,"f-",l,"g+",l,"g-",l
1255 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1256 ffthet(n,m,l,i,j,k,iblock),&
1257 ffthet(m,n,l,i,j,k,iblock),&
1258 ggthet(n,m,l,i,j,k,iblock),&
1259 ggthet(m,n,l,i,j,k,iblock)
1270 !C here will be the apropriate recalibrating for D-aminoacid
1271 read (ithep,*,end=121,err=121) nthetyp
1272 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1273 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1274 do i=-nthetyp+1,nthetyp-1
1275 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1276 do j=0,nbend_kcc_Tb(i)
1277 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1281 write (iout,'(a)') &
1282 "Parameters of the valence-only potentials"
1283 do i=-nthetyp+1,nthetyp-1
1284 write (iout,'(2a)') "Type ",toronelet(i)
1285 do k=0,nbend_kcc_Tb(i)
1286 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1292 write (2,*) "Start reading THETA_PDB",ithep_pdb
1294 ! write (2,*) 'i=',i
1295 read (ithep_pdb,*,err=111,end=111) &
1296 a0thet(i),(athet(j,i,1,1),j=1,2),&
1297 (bthet(j,i,1,1),j=1,2)
1298 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1299 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1300 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1301 sigc0(i)=sigc0(i)**2
1304 athet(1,i,1,-1)=athet(1,i,1,1)
1305 athet(2,i,1,-1)=athet(2,i,1,1)
1306 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1307 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1308 athet(1,i,-1,1)=-athet(1,i,1,1)
1309 athet(2,i,-1,1)=-athet(2,i,1,1)
1310 bthet(1,i,-1,1)=bthet(1,i,1,1)
1311 bthet(2,i,-1,1)=bthet(2,i,1,1)
1314 a0thet(i)=a0thet(-i)
1315 athet(1,i,-1,-1)=athet(1,-i,1,1)
1316 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1317 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1318 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1319 athet(1,i,-1,1)=athet(1,-i,1,1)
1320 athet(2,i,-1,1)=-athet(2,-i,1,1)
1321 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1322 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1323 athet(1,i,1,-1)=-athet(1,-i,1,1)
1324 athet(2,i,1,-1)=athet(2,-i,1,1)
1325 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1326 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1327 theta0(i)=theta0(-i)
1331 polthet(j,i)=polthet(j,-i)
1334 gthet(j,i)=gthet(j,-i)
1337 write (2,*) "End reading THETA_PDB"
1341 !--------------- Reading theta parameters for nucleic acid-------
1342 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1343 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1344 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1345 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1346 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1347 nthetyp_nucl+1,nthetyp_nucl+1))
1348 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1349 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1350 nthetyp_nucl+1,nthetyp_nucl+1))
1351 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1352 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1353 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1354 nthetyp_nucl+1,nthetyp_nucl+1))
1355 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1356 nthetyp_nucl+1,nthetyp_nucl+1))
1357 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1358 nthetyp_nucl+1,nthetyp_nucl+1))
1359 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1360 nthetyp_nucl+1,nthetyp_nucl+1))
1361 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1362 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1363 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1364 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1365 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1366 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1368 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1369 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1371 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1373 aa0thet_nucl(:,:,:)=0.0d0
1374 aathet_nucl(:,:,:,:)=0.0d0
1375 bbthet_nucl(:,:,:,:,:)=0.0d0
1376 ccthet_nucl(:,:,:,:,:)=0.0d0
1377 ddthet_nucl(:,:,:,:,:)=0.0d0
1378 eethet_nucl(:,:,:,:,:)=0.0d0
1379 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1380 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1385 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1386 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1387 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1388 read (ithep_nucl,*,end=111,err=111) &
1389 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1390 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1391 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1392 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1393 read (ithep_nucl,*,end=111,err=111) &
1394 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1395 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1396 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1401 !-------------------------------------------
1402 allocate(nlob(ntyp1)) !(ntyp1)
1403 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1404 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1405 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1412 gaussc(:,:,:,:)=0.0D0
1416 ! Read the parameters of the probability distribution/energy expression
1417 ! of the side chains.
1420 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1424 dsc_inv(i)=1.0D0/dsc(i)
1435 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1436 ((blower(k,l,1),l=1,k),k=1,3)
1437 censc(1,1,-i)=censc(1,1,i)
1438 censc(2,1,-i)=censc(2,1,i)
1439 censc(3,1,-i)=-censc(3,1,i)
1441 read (irotam,*,end=112,err=112) bsc(j,i)
1442 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1443 ((blower(k,l,j),l=1,k),k=1,3)
1444 censc(1,j,-i)=censc(1,j,i)
1445 censc(2,j,-i)=censc(2,j,i)
1446 censc(3,j,-i)=-censc(3,j,i)
1447 ! BSC is amplitude of Gaussian
1454 akl=akl+blower(k,m,j)*blower(l,m,j)
1458 if (((k.eq.3).and.(l.ne.3)) &
1459 .or.((l.eq.3).and.(k.ne.3))) then
1460 gaussc(k,l,j,-i)=-akl
1461 gaussc(l,k,j,-i)=-akl
1463 gaussc(k,l,j,-i)=akl
1464 gaussc(l,k,j,-i)=akl
1473 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1476 if (nlobi.gt.0) then
1478 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1479 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1480 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1481 'log h',(bsc(j,i),j=1,nlobi)
1482 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1483 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1485 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1486 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1489 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1490 write (iout,'(a,f10.4,4(16x,f10.4))') &
1491 'Center ',(bsc(j,i),j=1,nlobi)
1492 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1501 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1502 ! added by Urszula Kozlowska 07/11/2007
1504 !el Maximum number of SC local term fitting function coefficiants
1505 !el integer,parameter :: maxsccoef=65
1507 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1510 read (irotam,*,end=112,err=112)
1512 read (irotam,*,end=112,err=112)
1515 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1519 !---------reading nucleic acid parameters for rotamers-------------------
1520 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1521 do i=1,ntyp_molec(2)
1522 read (irotam_nucl,*,end=112,err=112)
1524 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1530 write (iout,*) "Base rotamer parameters"
1531 do i=1,ntyp_molec(2)
1532 write (iout,'(a)') restyp(i,2)
1533 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1538 ! Read the parameters of the probability distribution/energy expression
1539 ! of the side chains.
1541 write (2,*) "Start reading ROTAM_PDB"
1543 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1547 dsc_inv(i)=1.0D0/dsc(i)
1558 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1559 ((blower(k,l,1),l=1,k),k=1,3)
1561 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1562 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1563 ((blower(k,l,j),l=1,k),k=1,3)
1570 akl=akl+blower(k,m,j)*blower(l,m,j)
1580 write (2,*) "End reading ROTAM_PDB"
1586 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1587 !C interaction energy of the Gly, Ala, and Pro prototypes.
1589 read (ifourier,*) nloctyp
1590 SPLIT_FOURIERTOR = nloctyp.lt.0
1591 nloctyp = iabs(nloctyp)
1592 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1593 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1594 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1595 !C allocate(ctilde(2,2,nres))
1596 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1597 !C allocate(gtb1(2,nres))
1598 !C allocate(gtb2(2,nres))
1599 !C allocate(cc(2,2,nres))
1600 !C allocate(dd(2,2,nres))
1601 !C allocate(ee(2,2,nres))
1602 !C allocate(gtcc(2,2,nres))
1603 !C allocate(gtdd(2,2,nres))
1604 !C allocate(gtee(2,2,nres))
1607 allocate(itype2loc(-ntyp1:ntyp1))
1608 allocate(iloctyp(-nloctyp:nloctyp))
1609 allocate(bnew1(3,2,-nloctyp:nloctyp))
1610 allocate(bnew2(3,2,-nloctyp:nloctyp))
1611 allocate(ccnew(3,2,-nloctyp:nloctyp))
1612 allocate(ddnew(3,2,-nloctyp:nloctyp))
1613 allocate(e0new(3,-nloctyp:nloctyp))
1614 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1615 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1616 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1617 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1618 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1619 allocate(e0newtor(3,-nloctyp:nloctyp))
1620 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1633 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1634 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1635 itype2loc(ntyp1)=nloctyp
1636 iloctyp(nloctyp)=ntyp1
1638 itype2loc(-i)=-itype2loc(i)
1641 allocate(iloctyp(-nloctyp:nloctyp))
1642 allocate(itype2loc(-ntyp1:ntyp1))
1649 iloctyp(-i)=-iloctyp(i)
1651 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1652 !c write (iout,*) "nloctyp",nloctyp,
1653 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1654 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1655 !c write (iout,*) "nloctyp",nloctyp,
1656 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1659 !c write (iout,*) "NEWCORR",i
1660 read (ifourier,*,end=115,err=115)
1663 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1666 !c write (iout,*) "NEWCORR BNEW1"
1667 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1670 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1673 !c write (iout,*) "NEWCORR BNEW2"
1674 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1676 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1677 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1679 !c write (iout,*) "NEWCORR CCNEW"
1680 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1682 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1683 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1685 !c write (iout,*) "NEWCORR DDNEW"
1686 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1690 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1694 !c write (iout,*) "NEWCORR EENEW1"
1695 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1697 read (ifourier,*,end=115,err=115) e0new(ii,i)
1699 !c write (iout,*) (e0new(ii,i),ii=1,3)
1701 !c write (iout,*) "NEWCORR EENEW"
1704 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1705 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1706 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1707 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1712 bnew1(ii,1,-i)= bnew1(ii,1,i)
1713 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1714 bnew2(ii,1,-i)= bnew2(ii,1,i)
1715 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1718 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1719 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1720 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1721 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1722 ccnew(ii,1,-i)=ccnew(ii,1,i)
1723 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1724 ddnew(ii,1,-i)=ddnew(ii,1,i)
1725 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1727 e0new(1,-i)= e0new(1,i)
1728 e0new(2,-i)=-e0new(2,i)
1729 e0new(3,-i)=-e0new(3,i)
1731 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1732 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1733 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1734 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1738 write (iout,'(a)') "Coefficients of the multibody terms"
1739 do i=-nloctyp+1,nloctyp-1
1740 write (iout,*) "Type: ",onelet(iloctyp(i))
1741 write (iout,*) "Coefficients of the expansion of B1"
1743 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1745 write (iout,*) "Coefficients of the expansion of B2"
1747 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1749 write (iout,*) "Coefficients of the expansion of C"
1750 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1751 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1752 write (iout,*) "Coefficients of the expansion of D"
1753 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1754 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1755 write (iout,*) "Coefficients of the expansion of E"
1756 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1759 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1764 IF (SPLIT_FOURIERTOR) THEN
1766 !c write (iout,*) "NEWCORR TOR",i
1767 read (ifourier,*,end=115,err=115)
1770 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1773 !c write (iout,*) "NEWCORR BNEW1 TOR"
1774 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1777 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1780 !c write (iout,*) "NEWCORR BNEW2 TOR"
1781 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1783 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1784 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1786 !c write (iout,*) "NEWCORR CCNEW TOR"
1787 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1789 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1790 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1792 !c write (iout,*) "NEWCORR DDNEW TOR"
1793 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1797 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1801 !c write (iout,*) "NEWCORR EENEW1 TOR"
1802 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1804 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1806 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1808 !c write (iout,*) "NEWCORR EENEW TOR"
1811 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1812 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1813 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1814 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1819 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1820 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1821 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1822 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1825 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1826 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1827 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1828 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1830 e0newtor(1,-i)= e0newtor(1,i)
1831 e0newtor(2,-i)=-e0newtor(2,i)
1832 e0newtor(3,-i)=-e0newtor(3,i)
1834 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1835 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1836 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1837 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1841 write (iout,'(a)') &
1842 "Single-body coefficients of the torsional potentials"
1843 do i=-nloctyp+1,nloctyp-1
1844 write (iout,*) "Type: ",onelet(iloctyp(i))
1845 write (iout,*) "Coefficients of the expansion of B1tor"
1847 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1848 j,(bnew1tor(k,j,i),k=1,3)
1850 write (iout,*) "Coefficients of the expansion of B2tor"
1852 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1853 j,(bnew2tor(k,j,i),k=1,3)
1855 write (iout,*) "Coefficients of the expansion of Ctor"
1856 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1857 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1858 write (iout,*) "Coefficients of the expansion of Dtor"
1859 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1860 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1861 write (iout,*) "Coefficients of the expansion of Etor"
1862 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1865 write (iout,'(1hE,2i1,2f10.5)') &
1866 j,k,(eenewtor(l,j,k,i),l=1,2)
1872 do i=-nloctyp+1,nloctyp-1
1875 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1880 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1884 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1885 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1886 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1887 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1890 ENDIF !SPLIT_FOURIER_TOR
1892 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1893 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1894 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1895 allocate(b(13,-nloctyp-1:nloctyp+1))
1897 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1899 read (ifourier,*,end=115,err=115)
1900 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1902 write (iout,*) 'Type ',onelet(iloctyp(i))
1903 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1917 !c B1tilde(1,i) = b(3)
1918 !c! B1tilde(2,i) =-b(5)
1919 !c! B1tilde(1,-i) =-b(3)
1920 !c! B1tilde(2,-i) =b(5)
1921 !c! b1tilde(1,i)=0.0d0
1922 !c b1tilde(2,i)=0.0d0
1927 !cc B1tilde(1,i) = b(3,i)
1928 !cc B1tilde(2,i) =-b(5,i)
1929 !c B1tilde(1,-i) =-b(3,i)
1930 !c B1tilde(2,-i) =b(5,i)
1931 !cc b1tilde(1,i)=0.0d0
1932 !cc b1tilde(2,i)=0.0d0
1933 !cc B2(1,i) = b(2,i)
1934 !cc B2(2,i) = b(4,i)
1936 !c B2(2,-i) =-b(4,i)
1940 CCold(1,1,i)= b(7,i)
1941 CCold(2,2,i)=-b(7,i)
1942 CCold(2,1,i)= b(9,i)
1943 CCold(1,2,i)= b(9,i)
1944 CCold(1,1,-i)= b(7,i)
1945 CCold(2,2,-i)=-b(7,i)
1946 CCold(2,1,-i)=-b(9,i)
1947 CCold(1,2,-i)=-b(9,i)
1952 !c Ctilde(1,1,i)= CCold(1,1,i)
1953 !c Ctilde(1,2,i)= CCold(1,2,i)
1954 !c Ctilde(2,1,i)=-CCold(2,1,i)
1955 !c Ctilde(2,2,i)=-CCold(2,2,i)
1960 !c Ctilde(1,1,i)= CCold(1,1,i)
1961 !c Ctilde(1,2,i)= CCold(1,2,i)
1962 !c Ctilde(2,1,i)=-CCold(2,1,i)
1963 !c Ctilde(2,2,i)=-CCold(2,2,i)
1965 !c Ctilde(1,1,i)=0.0d0
1966 !c Ctilde(1,2,i)=0.0d0
1967 !c Ctilde(2,1,i)=0.0d0
1968 !c Ctilde(2,2,i)=0.0d0
1969 DDold(1,1,i)= b(6,i)
1970 DDold(2,2,i)=-b(6,i)
1971 DDold(2,1,i)= b(8,i)
1972 DDold(1,2,i)= b(8,i)
1973 DDold(1,1,-i)= b(6,i)
1974 DDold(2,2,-i)=-b(6,i)
1975 DDold(2,1,-i)=-b(8,i)
1976 DDold(1,2,-i)=-b(8,i)
1981 !c Dtilde(1,1,i)= DD(1,1,i)
1982 !c Dtilde(1,2,i)= DD(1,2,i)
1983 !c Dtilde(2,1,i)=-DD(2,1,i)
1984 !c Dtilde(2,2,i)=-DD(2,2,i)
1986 !c Dtilde(1,1,i)=0.0d0
1987 !c Dtilde(1,2,i)=0.0d0
1988 !c Dtilde(2,1,i)=0.0d0
1989 !c Dtilde(2,2,i)=0.0d0
1990 EEold(1,1,i)= b(10,i)+b(11,i)
1991 EEold(2,2,i)=-b(10,i)+b(11,i)
1992 EEold(2,1,i)= b(12,i)-b(13,i)
1993 EEold(1,2,i)= b(12,i)+b(13,i)
1994 EEold(1,1,-i)= b(10,i)+b(11,i)
1995 EEold(2,2,-i)=-b(10,i)+b(11,i)
1996 EEold(2,1,-i)=-b(12,i)+b(13,i)
1997 EEold(1,2,-i)=-b(12,i)-b(13,i)
1998 write(iout,*) "TU DOCHODZE"
2004 !c ee(2,1,i)=ee(1,2,i)
2009 "Coefficients of the cumulants (independent of valence angles)"
2010 do i=-nloctyp+1,nloctyp-1
2011 write (iout,*) 'Type ',onelet(iloctyp(i))
2013 write(iout,'(2f10.5)') B(3,i),B(5,i)
2015 write(iout,'(2f10.5)') B(2,i),B(4,i)
2018 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
2022 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
2026 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
2035 ! Read torsional parameters in old format
2037 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
2039 read (itorp,*,end=113,err=113) ntortyp,nterm_old
2040 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
2041 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2043 !el from energy module--------
2044 allocate(v1(nterm_old,ntortyp,ntortyp))
2045 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2046 !el---------------------------
2051 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2057 write (iout,'(/a/)') 'Torsional constants:'
2060 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2061 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2067 ! Read torsional parameters
2069 IF (TOR_MODE.eq.0) THEN
2070 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2071 read (itorp,*,end=113,err=113) ntortyp
2072 !el from energy module---------
2073 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2074 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2076 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2077 allocate(vlor2(maxlor,ntortyp,ntortyp))
2078 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2079 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2081 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2082 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2083 !el---------------------------
2086 !el---------------------------
2088 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2090 itortyp(i)=-itortyp(-i)
2092 itortyp(ntyp1)=ntortyp
2093 itortyp(-ntyp1)=-ntortyp
2095 write (iout,*) 'ntortyp',ntortyp
2097 do j=-ntortyp+1,ntortyp-1
2098 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2100 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2101 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2104 do k=1,nterm(i,j,iblock)
2105 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2107 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2108 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2109 v0ij=v0ij+si*v1(k,i,j,iblock)
2111 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2112 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2113 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2115 do k=1,nlor(i,j,iblock)
2116 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2117 vlor2(k,i,j),vlor3(k,i,j)
2118 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2121 v0(-i,-j,iblock)=v0ij
2127 write (iout,'(/a/)') 'Torsional constants:'
2129 do i=-ntortyp,ntortyp
2130 do j=-ntortyp,ntortyp
2131 write (iout,*) 'ityp',i,' jtyp',j
2132 write (iout,*) 'Fourier constants'
2133 do k=1,nterm(i,j,iblock)
2134 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2137 write (iout,*) 'Lorenz constants'
2138 do k=1,nlor(i,j,iblock)
2139 write (iout,'(3(1pe15.5))') &
2140 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2146 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2148 ! 6/23/01 Read parameters for double torsionals
2150 !el from energy module------------
2151 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2152 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2153 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2154 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2155 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2156 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2157 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2158 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2159 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2160 !---------------------------------
2164 do j=-ntortyp+1,ntortyp-1
2165 do k=-ntortyp+1,ntortyp-1
2166 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2167 ! write (iout,*) "OK onelett",
2170 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2171 .or. t3.ne.toronelet(k)) then
2172 write (iout,*) "Error in double torsional parameter file",&
2175 call MPI_Finalize(Ierror)
2177 stop "Error in double torsional parameter file"
2179 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2180 ntermd_2(i,j,k,iblock)
2181 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2182 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2183 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2184 ntermd_1(i,j,k,iblock))
2185 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2186 ntermd_1(i,j,k,iblock))
2187 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2188 ntermd_1(i,j,k,iblock))
2189 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2190 ntermd_1(i,j,k,iblock))
2191 ! Martix of D parameters for one dimesional foureir series
2192 do l=1,ntermd_1(i,j,k,iblock)
2193 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2194 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2195 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2196 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2197 ! write(iout,*) "whcodze" ,
2198 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2200 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2201 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2202 v2s(m,l,i,j,k,iblock),&
2203 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2204 ! Martix of D parameters for two dimesional fourier series
2205 do l=1,ntermd_2(i,j,k,iblock)
2207 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2208 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2209 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2210 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2219 write (iout,*) 'Constants for double torsionals'
2222 do j=-ntortyp+1,ntortyp-1
2223 do k=-ntortyp+1,ntortyp-1
2224 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2225 ' nsingle',ntermd_1(i,j,k,iblock),&
2226 ' ndouble',ntermd_2(i,j,k,iblock)
2228 write (iout,*) 'Single angles:'
2229 do l=1,ntermd_1(i,j,k,iblock)
2230 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2231 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2232 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2233 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2236 write (iout,*) 'Pairs of angles:'
2237 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2238 do l=1,ntermd_2(i,j,k,iblock)
2239 write (iout,'(i5,20f10.5)') &
2240 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2243 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2244 do l=1,ntermd_2(i,j,k,iblock)
2245 write (iout,'(i5,20f10.5)') &
2246 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2247 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2257 itype2loc(i)=itortyp(i)
2261 ELSE IF (TOR_MODE.eq.1) THEN
2263 !C read valence-torsional parameters
2264 read (itorp,*,end=121,err=121) ntortyp
2266 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2268 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2269 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2272 itype2loc(i)=itortyp(i)
2276 itortyp(i)=-itortyp(-i)
2278 do i=-ntortyp+1,ntortyp-1
2279 do j=-ntortyp+1,ntortyp-1
2280 !C first we read the cos and sin gamma parameters
2281 read (itorp,'(13x,a)',end=121,err=121) string
2282 write (iout,*) i,j,string
2283 read (itorp,*,end=121,err=121) &
2284 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2285 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2286 do k=1,nterm_kcc(j,i)
2287 do l=1,nterm_kcc_Tb(j,i)
2288 do ll=1,nterm_kcc_Tb(j,i)
2289 read (itorp,*,end=121,err=121) ii,jj,kk, &
2290 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2298 !c AL 4/8/16: Calculate coefficients from one-body parameters
2300 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2301 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2302 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2303 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2304 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2307 print *,i,itortyp(i)
2308 itortyp(i)=itype2loc(i)
2311 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2313 do i=-ntortyp+1,ntortyp-1
2314 do j=-ntortyp+1,ntortyp-1
2317 do k=1,nterm_kcc_Tb(j,i)
2318 do l=1,nterm_kcc_Tb(j,i)
2319 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2320 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2321 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2322 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2325 do k=1,nterm_kcc_Tb(j,i)
2326 do l=1,nterm_kcc_Tb(j,i)
2328 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2329 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2330 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2331 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2333 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2334 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2335 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2336 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2340 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2344 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2349 if (tor_mode.gt.0 .and. lprint) then
2350 !c Print valence-torsional parameters
2351 write (iout,'(a)') &
2352 "Parameters of the valence-torsional potentials"
2353 do i=-ntortyp+1,ntortyp-1
2354 do j=-ntortyp+1,ntortyp-1
2355 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2356 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2357 do k=1,nterm_kcc(j,i)
2358 do l=1,nterm_kcc_Tb(j,i)
2359 do ll=1,nterm_kcc_Tb(j,i)
2360 write (iout,'(3i5,2f15.4)')&
2361 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2369 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2370 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2371 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2372 !el from energy module---------
2373 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2374 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2376 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2377 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2378 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2379 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2381 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2382 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2383 !el---------------------------
2386 !el--------------------
2387 read (itorp_nucl,*,end=113,err=113) &
2388 (itortyp_nucl(i),i=1,ntyp_molec(2))
2389 ! print *,itortyp_nucl(:)
2390 !c write (iout,*) 'ntortyp',ntortyp
2393 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2394 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2397 do k=1,nterm_nucl(i,j)
2398 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2399 v0ij=v0ij+si*v1_nucl(k,i,j)
2402 do k=1,nlor_nucl(i,j)
2403 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2404 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2405 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2411 ! Read of Side-chain backbone correlation parameters
2412 ! Modified 11 May 2012 by Adasko
2415 read (isccor,*,end=119,err=119) nsccortyp
2417 !el from module energy-------------
2418 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2419 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2420 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2421 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2422 !-----------------------------------
2424 !el from module energy-------------
2425 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2427 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2429 isccortyp(i)=-isccortyp(-i)
2431 iscprol=isccortyp(20)
2432 ! write (iout,*) 'ntortyp',ntortyp
2434 !c maxinter is maximum interaction sites
2435 !el from module energy---------
2436 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2437 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2438 -nsccortyp:nsccortyp))
2439 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2440 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2441 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2442 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2443 !-----------------------------------
2447 read (isccor,*,end=119,err=119) &
2448 nterm_sccor(i,j),nlor_sccor(i,j)
2454 nterm_sccor(-i,j)=nterm_sccor(i,j)
2455 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2456 nterm_sccor(i,-j)=nterm_sccor(i,j)
2457 do k=1,nterm_sccor(i,j)
2458 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2460 if (j.eq.iscprol) then
2461 if (i.eq.isccortyp(10)) then
2462 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2463 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2465 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2466 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2467 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2468 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2469 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2470 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2471 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2472 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2475 if (i.eq.isccortyp(10)) then
2476 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2477 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2479 if (j.eq.isccortyp(10)) then
2480 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2481 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2483 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2484 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2485 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2486 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2487 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2488 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2492 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2493 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2494 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2495 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2498 do k=1,nlor_sccor(i,j)
2499 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2500 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2501 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2502 (1+vlor3sccor(k,i,j)**2)
2504 v0sccor(l,i,j)=v0ijsccor
2505 v0sccor(l,-i,j)=v0ijsccor1
2506 v0sccor(l,i,-j)=v0ijsccor2
2507 v0sccor(l,-i,-j)=v0ijsccor3
2513 !el from module energy-------------
2514 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2516 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2517 ! write (iout,*) 'ntortyp',ntortyp
2519 !c maxinter is maximum interaction sites
2520 !el from module energy---------
2521 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2522 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2523 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2524 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2525 !-----------------------------------
2529 read (isccor,*,end=119,err=119) &
2530 nterm_sccor(i,j),nlor_sccor(i,j)
2534 do k=1,nterm_sccor(i,j)
2535 print *,"test",i,j,k,l
2536 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2538 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2541 do k=1,nlor_sccor(i,j)
2542 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2543 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2544 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2545 (1+vlor3sccor(k,i,j)**2)
2547 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2556 write (iout,'(/a/)') 'Torsional constants:'
2559 write (iout,*) 'ityp',i,' jtyp',j
2560 write (iout,*) 'Fourier constants'
2561 do k=1,nterm_sccor(i,j)
2562 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2564 write (iout,*) 'Lorenz constants'
2565 do k=1,nlor_sccor(i,j)
2566 write (iout,'(3(1pe15.5))') &
2567 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2574 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2575 ! interaction energy of the Gly, Ala, and Pro prototypes.
2578 ! Read electrostatic-interaction parameters
2583 write (iout,'(/a)') 'Electrostatic interaction constants:'
2584 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2585 'IT','JT','APP','BPP','AEL6','AEL3'
2587 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2588 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2589 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2590 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2595 app (i,j)=epp(i,j)*rri*rri
2596 bpp (i,j)=-2.0D0*epp(i,j)*rri
2597 ael6(i,j)=elpp6(i,j)*4.2D0**6
2598 ael3(i,j)=elpp3(i,j)*4.2D0**3
2600 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2606 ! Read side-chain interaction parameters.
2608 !el from module energy - COMMON.INTERACT-------
2609 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2610 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2611 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2613 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2614 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2615 allocate(epslip(ntyp,ntyp))
2623 !--------------------------------
2625 read (isidep,*,end=117,err=117) ipot,expon
2626 if (ipot.lt.1 .or. ipot.gt.5) then
2627 write (iout,'(2a)') 'Error while reading SC interaction',&
2628 'potential file - unknown potential type.'
2630 call MPI_Finalize(Ierror)
2636 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2637 ', exponents are ',expon,2*expon
2638 ! goto (10,20,30,30,40) ipot
2640 !----------------------- LJ potential ---------------------------------
2642 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2643 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2644 (sigma0(i),i=1,ntyp)
2646 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2647 write (iout,'(a/)') 'The epsilon array:'
2648 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2649 write (iout,'(/a)') 'One-body parameters:'
2650 write (iout,'(a,4x,a)') 'residue','sigma'
2651 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2654 !----------------------- LJK potential --------------------------------
2656 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2657 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2658 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2660 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2661 write (iout,'(a/)') 'The epsilon array:'
2662 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2663 write (iout,'(/a)') 'One-body parameters:'
2664 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2665 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2669 !---------------------- GB or BP potential -----------------------------
2672 ! print *,"I AM in SCELE",scelemode
2673 if (scelemode.eq.0) then
2675 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2677 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2678 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2679 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2680 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2682 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2685 ! For the GB potential convert sigma'**2 into chi'
2688 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2692 write (iout,'(/a/)') 'Parameters of the BP potential:'
2693 write (iout,'(a/)') 'The epsilon array:'
2694 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2695 write (iout,'(/a)') 'One-body parameters:'
2696 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2698 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2699 chip(i),alp(i),i=1,ntyp)
2702 ! print *,ntyp,"NTYP"
2703 allocate(icharge(ntyp1))
2704 ! print *,ntyp,icharge(i)
2706 read (isidep,*) (icharge(i),i=1,ntyp)
2707 print *,ntyp,icharge(i)
2708 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2709 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2710 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2711 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2712 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2713 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2714 allocate(epsintab(ntyp,ntyp))
2715 allocate(dtail(2,ntyp,ntyp))
2716 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2717 allocate(wqdip(2,ntyp,ntyp))
2718 allocate(wstate(4,ntyp,ntyp))
2719 allocate(dhead(2,2,ntyp,ntyp))
2720 allocate(nstate(ntyp,ntyp))
2721 allocate(debaykap(ntyp,ntyp))
2723 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2724 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2728 ! write (*,*) "Im in ALAB", i, " ", j
2730 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2731 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2732 chis(i,j),chis(j,i), & !2 w tej linii
2733 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2734 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2735 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2736 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2737 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2738 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2739 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2740 IF ((LaTeX).and.(i.gt.24)) then
2741 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2742 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2743 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2744 chis(i,j),chis(j,i) !2 w tej linii
2746 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2750 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2754 IF ((LaTeX).and.(i.gt.24)) then
2755 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2756 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2757 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2758 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2759 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2760 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2761 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2768 sigma(i,j) = sigma(j,i)
2769 sigmap1(i,j)=sigmap1(j,i)
2770 sigmap2(i,j)=sigmap2(j,i)
2771 sigiso1(i,j)=sigiso1(j,i)
2772 sigiso2(i,j)=sigiso2(j,i)
2773 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2774 nstate(i,j) = nstate(j,i)
2775 dtail(1,i,j) = dtail(1,j,i)
2776 dtail(2,i,j) = dtail(2,j,i)
2778 alphasur(k,i,j) = alphasur(k,j,i)
2779 wstate(k,i,j) = wstate(k,j,i)
2780 alphiso(k,i,j) = alphiso(k,j,i)
2783 dhead(2,1,i,j) = dhead(1,1,j,i)
2784 dhead(2,2,i,j) = dhead(1,2,j,i)
2785 dhead(1,1,i,j) = dhead(2,1,j,i)
2786 dhead(1,2,i,j) = dhead(2,2,j,i)
2788 epshead(i,j) = epshead(j,i)
2789 sig0head(i,j) = sig0head(j,i)
2792 wqdip(k,i,j) = wqdip(k,j,i)
2795 wquad(i,j) = wquad(j,i)
2796 epsintab(i,j) = epsintab(j,i)
2797 debaykap(i,j)=debaykap(j,i)
2798 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2805 !--------------------- GBV potential -----------------------------------
2807 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2808 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2809 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2810 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2812 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2813 write (iout,'(a/)') 'The epsilon array:'
2814 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2815 write (iout,'(/a)') 'One-body parameters:'
2816 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2817 's||/s_|_^2',' chip ',' alph '
2818 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2819 sigii(i),chip(i),alp(i),i=1,ntyp)
2822 write(iout,*)"Wrong ipot"
2828 !-----------------------------------------------------------------------
2829 ! Calculate the "working" parameters of SC interactions.
2831 !el from module energy - COMMON.INTERACT-------
2832 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2833 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2834 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2835 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2836 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2837 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2839 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2845 if (scelemode.eq.0) then
2854 sc_aa_tube_par(:)=0.0d0
2855 sc_bb_tube_par(:)=0.0d0
2857 !--------------------------------
2862 epslip(i,j)=epslip(j,i)
2865 if (scelemode.eq.0) then
2868 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2869 sigma(j,i)=sigma(i,j)
2870 rs0(i,j)=dwa16*sigma(i,j)
2875 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2876 'Working parameters of the SC interactions:',&
2877 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2882 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2884 ! print *,"SIGMA ZLA?",sigma(i,j)
2892 sigeps=dsign(1.0D0,epsij)
2894 aa_aq(i,j)=epsij*rrij*rrij
2895 ! print *,"ADASKO",epsij,rrij,expon
2896 bb_aq(i,j)=-sigeps*epsij*rrij
2897 aa_aq(j,i)=aa_aq(i,j)
2898 bb_aq(j,i)=bb_aq(i,j)
2899 epsijlip=epslip(i,j)
2900 sigeps=dsign(1.0D0,epsijlip)
2901 epsijlip=dabs(epsijlip)
2902 aa_lip(i,j)=epsijlip*rrij*rrij
2903 bb_lip(i,j)=-sigeps*epsijlip*rrij
2904 aa_lip(j,i)=aa_lip(i,j)
2905 bb_lip(j,i)=bb_lip(i,j)
2906 !C write(iout,*) aa_lip
2907 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2908 sigt1sq=sigma0(i)**2
2909 sigt2sq=sigma0(j)**2
2912 ratsig1=sigt2sq/sigt1sq
2913 ratsig2=1.0D0/ratsig1
2914 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2915 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2916 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2920 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2921 sigmaii(i,j)=rsum_max
2922 sigmaii(j,i)=rsum_max
2924 ! sigmaii(i,j)=r0(i,j)
2925 ! sigmaii(j,i)=r0(i,j)
2927 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2928 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2929 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2930 augm(i,j)=epsij*r_augm**(2*expon)
2931 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2938 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2939 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2940 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2945 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2946 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2947 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2948 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2949 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2950 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2951 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2952 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2953 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2954 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2955 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2956 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2957 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2958 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2959 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2960 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2969 read (isidep_nucl,*) ipot_nucl
2970 ! print *,"TU?!",ipot_nucl
2971 if (ipot_nucl.eq.1) then
2972 do i=1,ntyp_molec(2)
2973 do j=i,ntyp_molec(2)
2974 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2975 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2979 do i=1,ntyp_molec(2)
2980 do j=i,ntyp_molec(2)
2981 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2982 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2983 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2987 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2988 do i=1,ntyp_molec(2)
2989 do j=i,ntyp_molec(2)
2990 rrij=sigma_nucl(i,j)
2994 epsij=4*eps_nucl(i,j)
2995 sigeps=dsign(1.0D0,epsij)
2997 aa_nucl(i,j)=epsij*rrij*rrij
2998 bb_nucl(i,j)=-sigeps*epsij*rrij
2999 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
3000 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
3001 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
3002 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
3003 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
3004 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
3007 aa_nucl(i,j)=aa_nucl(j,i)
3008 bb_nucl(i,j)=bb_nucl(j,i)
3009 ael3_nucl(i,j)=ael3_nucl(j,i)
3010 ael6_nucl(i,j)=ael6_nucl(j,i)
3011 ael63_nucl(i,j)=ael63_nucl(j,i)
3012 ael32_nucl(i,j)=ael32_nucl(j,i)
3013 elpp3_nucl(i,j)=elpp3_nucl(j,i)
3014 elpp6_nucl(i,j)=elpp6_nucl(j,i)
3015 elpp63_nucl(i,j)=elpp63_nucl(j,i)
3016 elpp32_nucl(i,j)=elpp32_nucl(j,i)
3017 eps_nucl(i,j)=eps_nucl(j,i)
3018 sigma_nucl(i,j)=sigma_nucl(j,i)
3019 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
3023 write(iout,*) "tube param"
3024 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
3025 ccavtubpep,dcavtubpep,tubetranenepep
3026 sigmapeptube=sigmapeptube**6
3027 sigeps=dsign(1.0D0,epspeptube)
3028 epspeptube=dabs(epspeptube)
3029 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
3030 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
3031 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
3033 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
3034 ccavtub(i),dcavtub(i),tubetranene(i)
3035 sigmasctube=sigmasctube**6
3036 sigeps=dsign(1.0D0,epssctube)
3037 epssctube=dabs(epssctube)
3038 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
3039 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
3040 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
3042 !-----------------READING SC BASE POTENTIALS-----------------------------
3043 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
3044 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
3045 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
3046 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3047 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3048 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3049 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3050 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3051 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3052 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3053 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3054 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3055 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3056 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3057 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3058 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3060 write (iout,*) "ESCBASEPARM"
3061 do i=1,ntyp_molec(1)
3062 do j=1,ntyp_molec(2) ! without U then we will take T for U
3063 ! write (*,*) "Im in ", i, " ", j
3064 read(isidep_scbase,*) &
3065 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3066 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3067 ! write(*,*) "eps",eps_scbase(i,j)
3068 read(isidep_scbase,*) &
3069 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3070 chis_scbase(i,j,1),chis_scbase(i,j,2)
3071 read(isidep_scbase,*) &
3072 dhead_scbasei(i,j), &
3073 dhead_scbasej(i,j), &
3074 rborn_scbasei(i,j),rborn_scbasej(i,j)
3075 read(isidep_scbase,*) &
3076 (wdipdip_scbase(k,i,j),k=1,3), &
3077 (wqdip_scbase(k,i,j),k=1,2)
3078 read(isidep_scbase,*) &
3079 alphapol_scbase(i,j), &
3080 epsintab_scbase(i,j)
3081 if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
3082 if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
3083 if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
3084 if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
3085 if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
3086 if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
3087 if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
3088 if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
3090 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3091 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3092 write(*,*) "eps",eps_scbase(i,j)
3094 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3095 chis_scbase(i,j,1),chis_scbase(i,j,2)
3097 dhead_scbasei(i,j), &
3098 dhead_scbasej(i,j), &
3099 rborn_scbasei(i,j),rborn_scbasej(i,j)
3101 (wdipdip_scbase(k,i,j),k=1,3), &
3102 (wqdip_scbase(k,i,j),k=1,2)
3104 alphapol_scbase(i,j), &
3105 epsintab_scbase(i,j)
3110 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3111 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3112 write(*,*) "eps",eps_scbase(i,j)
3114 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3115 chis_scbase(i,j,1),chis_scbase(i,j,2)
3117 dhead_scbasei(i,j), &
3118 dhead_scbasej(i,j), &
3119 rborn_scbasei(i,j),rborn_scbasej(i,j)
3121 (wdipdip_scbase(k,i,j),k=1,3), &
3122 (wqdip_scbase(k,i,j),k=1,2)
3124 alphapol_scbase(i,j), &
3125 epsintab_scbase(i,j)
3128 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3129 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3131 do i=1,ntyp_molec(1)
3132 do j=1,ntyp_molec(2)-1
3133 epsij=eps_scbase(i,j)
3134 rrij=sigma_scbase(i,j)
3139 sigeps=dsign(1.0D0,epsij)
3141 aa_scbase(i,j)=epsij*rrij*rrij
3142 bb_scbase(i,j)=-sigeps*epsij*rrij
3145 !-----------------READING PEP BASE POTENTIALS-------------------
3146 allocate(eps_pepbase(ntyp_molec(2)))
3147 allocate(sigma_pepbase(ntyp_molec(2)))
3148 allocate(chi_pepbase(ntyp_molec(2),2))
3149 allocate(chipp_pepbase(ntyp_molec(2),2))
3150 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3151 allocate(sigmap1_pepbase(ntyp_molec(2)))
3152 allocate(sigmap2_pepbase(ntyp_molec(2)))
3153 allocate(chis_pepbase(ntyp_molec(2),2))
3154 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3156 write (iout,*) "EPEPBASEPARM"
3158 do j=1,ntyp_molec(2) ! without U then we will take T for U
3159 write (*,*) "Im in ", i, " ", j
3160 read(isidep_pepbase,*) &
3161 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3162 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3163 if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
3164 if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
3165 if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
3166 if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
3167 if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
3168 if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
3169 if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
3170 if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
3172 write(*,*) "eps",eps_pepbase(j)
3173 read(isidep_pepbase,*) &
3174 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3175 chis_pepbase(j,1),chis_pepbase(j,2)
3176 read(isidep_pepbase,*) &
3177 (wdipdip_pepbase(k,j),k=1,3)
3178 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3179 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3181 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3182 chis_pepbase(j,1),chis_pepbase(j,2)
3184 (wdipdip_pepbase(k,j),k=1,3)
3188 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3189 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3191 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3192 chis_pepbase(j,1),chis_pepbase(j,2)
3194 (wdipdip_pepbase(k,j),k=1,3)
3196 allocate(aa_pepbase(ntyp_molec(2)))
3197 allocate(bb_pepbase(ntyp_molec(2)))
3199 do j=1,ntyp_molec(2)-1
3200 epsij=eps_pepbase(j)
3201 rrij=sigma_pepbase(j)
3206 sigeps=dsign(1.0D0,epsij)
3208 aa_pepbase(j)=epsij*rrij*rrij
3209 bb_pepbase(j)=-sigeps*epsij*rrij
3211 !--------------READING SC PHOSPHATE-------------------------------------
3212 allocate(eps_scpho(ntyp_molec(1)))
3213 allocate(sigma_scpho(ntyp_molec(1)))
3214 allocate(chi_scpho(ntyp_molec(1),2))
3215 allocate(chipp_scpho(ntyp_molec(1),2))
3216 allocate(alphasur_scpho(4,ntyp_molec(1)))
3217 allocate(sigmap1_scpho(ntyp_molec(1)))
3218 allocate(sigmap2_scpho(ntyp_molec(1)))
3219 allocate(chis_scpho(ntyp_molec(1),2))
3220 allocate(wqq_scpho(ntyp_molec(1)))
3221 allocate(wqdip_scpho(2,ntyp_molec(1)))
3222 allocate(alphapol_scpho(ntyp_molec(1)))
3223 allocate(epsintab_scpho(ntyp_molec(1)))
3224 allocate(dhead_scphoi(ntyp_molec(1)))
3225 allocate(rborn_scphoi(ntyp_molec(1)))
3226 allocate(rborn_scphoj(ntyp_molec(1)))
3227 allocate(alphi_scpho(ntyp_molec(1)))
3231 do j=1,ntyp_molec(1) ! without U then we will take T for U
3232 write (*,*) "Im in scpho ", i, " ", j
3233 read(isidep_scpho,*) &
3234 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3235 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3236 write(*,*) "eps",eps_scpho(j)
3237 read(isidep_scpho,*) &
3238 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3239 chis_scpho(j,1),chis_scpho(j,2)
3240 read(isidep_scpho,*) &
3241 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3242 read(isidep_scpho,*) &
3243 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3245 if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
3246 if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
3247 if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
3248 if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
3249 if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
3250 if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
3251 if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
3252 if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
3256 allocate(aa_scpho(ntyp_molec(1)))
3257 allocate(bb_scpho(ntyp_molec(1)))
3259 do j=1,ntyp_molec(1)
3266 sigeps=dsign(1.0D0,epsij)
3268 aa_scpho(j)=epsij*rrij*rrij
3269 bb_scpho(j)=-sigeps*epsij*rrij
3273 read(isidep_peppho,*) &
3274 eps_peppho,sigma_peppho
3275 read(isidep_peppho,*) &
3276 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3277 read(isidep_peppho,*) &
3278 (wqdip_peppho(k),k=1,2)
3286 sigeps=dsign(1.0D0,epsij)
3288 aa_peppho=epsij*rrij*rrij
3289 bb_peppho=-sigeps*epsij*rrij
3292 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3297 ! Define the SC-p interaction constants (hard-coded; old style)
3300 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3302 ! aad(i,1)=0.3D0*4.0D0**12
3303 ! Following line for constants currently implemented
3304 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3305 aad(i,1)=1.5D0*4.0D0**12
3306 ! aad(i,1)=0.17D0*5.6D0**12
3308 ! "Soft" SC-p repulsion
3310 ! Following line for constants currently implemented
3311 ! aad(i,1)=0.3D0*4.0D0**6
3312 ! "Hard" SC-p repulsion
3313 bad(i,1)=3.0D0*4.0D0**6
3314 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3323 ! 8/9/01 Read the SC-p interaction constants from file
3326 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3329 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3330 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3331 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3332 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3336 write (iout,*) "Parameters of SC-p interactions:"
3338 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3339 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3344 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3346 do i=1,ntyp_molec(2)
3347 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3349 do i=1,ntyp_molec(2)
3350 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3351 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3353 r0pp=1.12246204830937298142*5.16158
3359 ! Define the constants of the disulfide bridge
3363 ! Old arbitrary potential - commented out.
3368 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3369 ! energy surface of diethyl disulfide.
3370 ! A. Liwo and U. Kozlowska, 11/24/03
3386 allocate(ichargelipid(ntyp_molec(4)))
3387 allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3388 allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3389 allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
3390 allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
3391 allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
3392 kjtokcal=0.2390057361
3394 !HERE THE MASS of MARTINI
3395 write(*,*) "before MARTINI PARAM"
3396 do i=1,ntyp_molec(4)
3402 !relative dielectric constant = 15 for implicit screening
3403 k_coulomb_lip=332.0d0/15.0d0
3404 !kbond = 1250 kJ/(mol*nm*2)
3405 kbondlip=1250.0d0*kjtokcal/100.0d0
3407 lip_angle_force=0.0d0
3408 if (DRY_MARTINI.gt.0) then
3409 lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
3410 lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
3411 lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
3412 lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
3413 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3414 lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
3415 lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
3416 lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
3418 lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
3419 lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
3420 lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
3421 lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
3422 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3423 lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
3424 lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
3425 lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
3427 lip_angle_angle=0.0d0
3428 lip_angle_angle(3,12,12)=120.0/krad
3429 lip_angle_angle(3,12,18)=180.0/krad
3430 lip_angle_angle(3,18,16)=180.0/krad
3431 lip_angle_angle(12,18,16)=180.0/krad
3432 lip_angle_angle(18,16,18)=120.0/krad
3433 lip_angle_angle(16,18,18)=180.0/krad
3434 lip_angle_angle(12,18,18)=180.0/krad
3435 lip_angle_angle(18,18,18)=180.0/krad
3436 read(ilipbond,*) temp1
3438 read(ilipbond,*) temp1, lip_bond(i,1), &
3439 lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
3440 lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
3441 lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
3442 lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
3445 lip_bond(i,j)=lip_bond(i,j)*10
3449 read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
3450 read(ilipnonbond,*) temp1
3452 read(ilipnonbond,*) temp1, lip_eps(i,1), &
3453 lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
3454 lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
3455 lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
3456 lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
3458 ! write(*,*) i, lip_eps(i,18)
3460 lip_eps(i,j)=lip_eps(i,j)*kjtokcal
3463 read(ilipnonbond,*) temp1
3465 read(ilipnonbond,*) temp1,lip_sig(i,1), &
3466 lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
3467 lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
3468 lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
3469 lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
3472 lip_sig(i,j)=lip_sig(i,j)*10.0
3475 write(*,*) "after MARTINI PARAM"
3479 allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
3480 allocate(alphapolcat2(ntyp,ntyp))
3481 allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
3482 allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
3483 allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
3484 allocate(epsintabcat(ntyp,ntyp))
3485 allocate(dtailcat(2,ntyp,ntyp))
3486 allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
3487 allocate(wqdipcat(2,ntyp,ntyp))
3488 allocate(wstatecat(4,ntyp,ntyp))
3489 allocate(dheadcat(2,2,ntyp,ntyp))
3490 allocate(nstatecat(ntyp,ntyp))
3491 allocate(debaykapcat(ntyp,ntyp))
3493 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
3494 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
3495 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3496 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3497 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
3500 if (.not.allocated(ichargecat))&
3501 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
3502 write(*,*) "before ions",oldion
3505 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3506 if (oldion.eq.0) then
3507 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3508 allocate(icharge(1:ntyp1))
3509 read(iion,*) (icharge(i),i=1,ntyp)
3513 print *,ntyp_molec(5)
3514 do i=-ntyp_molec(5),ntyp_molec(5)
3515 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3516 print *,msc(i,5),restok(i,5)
3522 do j=1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
3524 ! do j=1,ntyp_molec(5)
3525 ! write (*,*) "Im in ALAB", i, " ", j
3527 epscat(i,j),sigmacat(i,j), &
3528 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3529 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
3531 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
3532 ! chiscat(i,j),chiscat(j,i), &
3533 chis1cat(i,j),chis2cat(i,j), &
3535 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
3536 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
3537 dtailcat(1,i,j),dtailcat(2,i,j), &
3538 epsheadcat(i,j),sig0headcat(i,j), &
3540 ! rborncat(i,j),rborncat(j,i),&
3541 rborn1cat(i,j),rborn2cat(i,j),&
3542 (wqdipcat(k,i,j),k=1,2), &
3543 alphapolcat(i,j),alphapolcat2(j,i), &
3544 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3546 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3547 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3549 ! write (iout,*) 'i= ', i, ' j= ', j
3550 ! write (iout,*) 'epsi0= ', epscat(1,j)
3551 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3552 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3553 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3554 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3555 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3556 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3557 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3558 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3559 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3560 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3561 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3562 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3563 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3564 ! write (iout,*) 'a1= ', rborncat(j,1)
3565 ! write (iout,*) 'a2= ', rborncat(1,j)
3566 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3567 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3568 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3569 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3573 ! If ((i.eq.1).and.(j.eq.27)) then
3574 ! write (iout,*) 'SEP'
3575 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3576 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3581 allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
3583 do j=1,ntyp_molec(5)
3587 sigeps=dsign(1.0D0,epsij)
3589 aa_aq_cat(i,j)=epsij*rrij*rrij
3590 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3595 do j=1,ntyp_molec(5)-1
3597 write (iout,*) 'i= ', i, ' j= ', j
3598 write (iout,*) 'epsi0= ', epscat(i,j)
3599 write (iout,*) 'sigma0= ', sigmacat(i,j)
3600 write (iout,*) 'chi1= ', chi1cat(i,j)
3601 write (iout,*) 'chi1= ', chi2cat(i,j)
3602 write (iout,*) 'chip1= ', chipp1cat(i,j)
3603 write (iout,*) 'chip2= ', chipp2cat(i,j)
3604 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3605 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3606 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3607 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3608 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3609 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3610 write (iout,*) 'chis1= ', chis1cat(i,j)
3611 write (iout,*) 'chis1= ', chis2cat(i,j)
3612 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3613 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3614 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3615 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3616 write (iout,*) 'a1= ', rborn1cat(i,j)
3617 write (iout,*) 'a2= ', rborn2cat(i,j)
3618 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3619 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3620 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3621 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3622 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3623 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3626 If ((i.eq.1).and.(j.eq.27)) then
3627 write (iout,*) 'SEP'
3628 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3629 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3636 ! read number of Zn2+
3637 ! here two denotes the Zn2+ and Cu2+
3638 write(iout,*) "before TRANPARM"
3639 allocate(aomicattr(0:3,2))
3640 allocate(athetacattran(0:6,5,2))
3641 allocate(agamacattran(3,5,2))
3642 allocate(acatshiftdsc(5,2))
3643 allocate(bcatshiftdsc(5,2))
3644 allocate(demorsecat(5,2))
3645 allocate(alphamorsecat(5,2))
3646 allocate(x0catleft(5,2))
3647 allocate(x0catright(5,2))
3648 allocate(x0cattrans(5,2))
3649 allocate(ntrantyp(2))
3650 do i=1,1 ! currently only Zn2+
3652 read(iiontran,*) ntrantyp(i)
3655 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3656 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
3657 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
3658 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3659 read(iiontran,*) (aomicattr(j,i),j=0,3)
3661 read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
3662 (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
3663 demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
3668 write (iout,'(/a)') "Disulfide bridge parameters:"
3669 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3670 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3671 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3672 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3675 if (shield_mode.gt.0) then
3676 pi=4.0D0*datan(1.0D0)
3677 !C VSolvSphere the volume of solving sphere
3679 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3680 !C there will be no distinction between proline peptide group and normal peptide
3681 !C group in case of shielding parameters
3682 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3683 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3684 write (iout,*) VSolvSphere,VSolvSphere_div
3685 !C long axis of side chain
3687 long_r_sidechain(i)=vbldsc0(1,i)
3688 ! if (scelemode.eq.0) then
3689 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3690 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3692 ! short_r_sidechain(i)=sigma(i,i)
3694 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3701 111 write (iout,*) "Error reading bending energy parameters."
3703 112 write (iout,*) "Error reading rotamer energy parameters."
3705 113 write (iout,*) "Error reading torsional energy parameters."
3707 114 write (iout,*) "Error reading double torsional energy parameters."
3709 115 write (iout,*) &
3710 "Error reading cumulant (multibody energy) parameters."
3712 116 write (iout,*) "Error reading electrostatic energy parameters."
3714 117 write (iout,*) "Error reading side chain interaction parameters."
3716 118 write (iout,*) "Error reading SCp interaction parameters."
3718 119 write (iout,*) "Error reading SCCOR parameters"
3720 121 write (iout,*) "Error in Czybyshev parameters"
3722 123 write(iout,*) "Error in transition metal parameters"
3725 call MPI_Finalize(Ierror)
3729 end subroutine parmread
3731 !-----------------------------------------------------------------------------
3733 !-----------------------------------------------------------------------------
3734 subroutine printmat(ldim,m,n,iout,key,a)
3737 character(len=3),dimension(n) :: key
3738 real(kind=8),dimension(ldim,n) :: a
3740 integer :: i,j,k,m,iout,nlim
3744 write (iout,1000) (key(k),k=i,nlim)
3746 1000 format (/5x,8(6x,a3))
3747 1020 format (/80(1h-)/)
3749 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3752 1010 format (a3,2x,8(f9.4))
3754 end subroutine printmat
3755 !-----------------------------------------------------------------------------
3757 !-----------------------------------------------------------------------------
3759 ! Read the PDB file and convert the peptide geometry into virtual-chain
3762 use energy_data, only: itype,istype
3766 ! use control, only: rescode,sugarcode
3767 ! implicit real*8 (a-h,o-z)
3768 ! include 'DIMENSIONS'
3769 ! include 'COMMON.LOCAL'
3770 ! include 'COMMON.VAR'
3771 ! include 'COMMON.CHAIN'
3772 ! include 'COMMON.INTERACT'
3773 ! include 'COMMON.IOUNITS'
3774 ! include 'COMMON.GEO'
3775 ! include 'COMMON.NAMES'
3776 ! include 'COMMON.CONTROL'
3777 ! include 'COMMON.DISTFIT'
3778 ! include 'COMMON.SETUP'
3779 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3781 logical :: lprn=.true.,fail
3782 real(kind=8),dimension(3) :: e1,e2,e3
3783 real(kind=8) :: dcj,efree_temp
3784 character(len=3) :: seq,res,res2
3785 character(len=5) :: atom
3786 character(len=80) :: card
3787 real(kind=8),dimension(3,40) :: sccor
3788 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3789 integer :: isugar,molecprev,firstion
3790 character*1 :: sugar
3792 real(kind=8),dimension(3) :: ccc
3794 integer,dimension(2,maxres/3) :: hfrag_alloc
3795 integer,dimension(4,maxres/3) :: bfrag_alloc
3796 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3797 real(kind=8),dimension(:,:), allocatable :: c_temporary
3798 integer,dimension(:,:) , allocatable :: itype_temporary
3799 integer,dimension(:),allocatable :: istype_temp
3806 ! write (2,*) "UNRES_PDB",unres_pdb
3826 !-----------------------------
3827 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3828 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3829 if(.not. allocated(istype)) allocate(istype(maxres))
3831 read (ipdbin,'(a80)',end=10) card
3832 write (iout,'(a)') card
3833 if (card(:5).eq.'HELIX') then
3836 read(card(22:25),*) hfrag(1,nhfrag)
3837 read(card(34:37),*) hfrag(2,nhfrag)
3839 if (card(:5).eq.'SHEET') then
3842 read(card(24:26),*) bfrag(1,nbfrag)
3843 read(card(35:37),*) bfrag(2,nbfrag)
3844 !rc----------------------------------------
3845 !rc to be corrected !!!
3846 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3847 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3848 !rc----------------------------------------
3850 if (card(:3).eq.'END') then
3852 else if (card(:3).eq.'TER') then
3857 itype(ires_old,molecule)=ntyp1_molec(molecule)
3858 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3859 nres_molec(molecule)=nres_molec(molecule)+2
3860 ! if (molecule.eq.4) ires=ires+2
3862 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3865 dc(j,ires)=sccor(j,iii)
3868 call sccenter(ires,iii,sccor)
3872 else if (card(:3).eq.'BRA') then
3876 itype(ires,molecule)=ntyp1_molec(molecule)-1
3877 nres_molec(molecule)=nres_molec(molecule)+1
3881 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3882 ! Fish out the ATOM cards.
3883 ! write(iout,*) 'card',card(1:20)
3884 ! print *,"ATU ",card(1:6), CARD(3:6)
3886 if (index(card(1:4),'ATOM').gt.0) then
3887 read (card(12:16),*) atom
3888 ! write (iout,*) "! ",atom," !",ires
3889 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3890 if (card(14:16).eq.'LIP') then
3895 nres_molec(molecule)=nres_molec(molecule)+1
3896 itype(ires,molecule)=ntyp1_molec(molecule)
3905 nres_molec(molecule)=nres_molec(molecule)+1
3906 read (card(18:20),'(a3)') res
3908 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3909 if (UNRES_PDB) then!
3910 itype(ires,molecule)=rescode(ires,res,0,molecule)
3912 itype(ires,molecule)=rescode_lip(res,ires)
3915 read (card(23:26),*) ires
3916 read (card(18:20),'(a3)') res
3917 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3918 ! & " ires_old",ires_old
3919 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3920 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3921 if (ires-ishift+ishift1.ne.ires_old) then
3922 ! Calculate the CM of the preceding residue.
3923 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3925 ! write (iout,*) "Calculating sidechain center iii",iii
3928 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3931 call sccenter(ires_old,iii,sccor)
3935 ! Start new residue.
3936 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3939 else if (ibeg.eq.1) then
3940 write (iout,*) "BEG ires",ires
3942 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3945 nres_molec(molecule)=nres_molec(molecule)+1
3947 ires=ires-ishift+ishift1
3949 ! write (iout,*) "ishift",ishift," ires",ires,&
3950 ! " ires_old",ires_old
3952 else if (ibeg.eq.2) then
3954 ishift=-ires_old+ires-1 !!!!!
3955 ishift1=ishift1-1 !!!!!
3956 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3957 ires=ires-ishift+ishift1
3958 ! print *,ires,ishift,ishift1
3962 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3963 ires=ires-ishift+ishift1
3966 ! print *,'atom',ires,atom
3967 if (res.eq.'ACE' .or. res.eq.'NHE') then
3970 if (atom.eq.'CA '.or.atom.eq.'N ') then
3972 itype(ires,molecule)=rescode(ires,res,0,molecule)
3974 ! nres_molec(molecule)=nres_molec(molecule)+1
3978 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3979 ! nres_molec(molecule)=nres_molec(molecule)+1
3980 read (card(19:19),'(a1)') sugar
3981 isugar=sugarcode(sugar,ires)
3982 ! if (ibeg.eq.1) then
3986 ! print *,"ires=",ires,istype(ires)
3992 ires=ires-ishift+ishift1
3994 ! write (iout,*) "ires_old",ires_old," ires",ires
3995 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3998 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3999 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
4000 res.eq.'NHE'.and.atom(:2).eq.'HN') then
4001 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4002 ! print *,ires,ishift,ishift1
4003 ! write (iout,*) "backbone ",atom
4005 write (iout,'(2i3,2x,a,3f8.3)') &
4006 ires,itype(ires,1),res,(c(j,ires),j=1,3)
4009 nres_molec(molecule)=nres_molec(molecule)+1
4011 sccor(j,iii)=c(j,ires)
4013 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
4014 atom.eq."C2'" .or. atom.eq."C3'" &
4015 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
4016 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
4017 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
4018 ! print *,ires,ishift,ishift1
4022 ! sccor(j,iii)=c(j,ires)
4025 c(j,ires)=c(j,ires)+ccc(j)/5.0
4027 print *,counter,molecule
4028 if (counter.eq.5) then
4030 nres_molec(molecule)=nres_molec(molecule)+1
4033 ! sccor(j,iii)=c(j,ires)
4037 ! print *, "ATOM",atom(1:3)
4038 else if (atom.eq."C5'") then
4039 read (card(19:19),'(a1)') sugar
4040 isugar=sugarcode(sugar,ires)
4045 ! print *,ires,istype(ires)
4049 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
4050 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4051 nres_molec(molecule)=nres_molec(molecule)+1
4052 print *,"nres_molec(molecule)",nres_molec(molecule),ires
4056 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4058 else if ((atom.eq."C1'").and.unres_pdb) then
4060 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4061 ! write (*,*) card(23:27),ires,itype(ires,1)
4062 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
4063 atom.ne.'N' .and. atom.ne.'C' .and. &
4064 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
4065 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
4066 .and. atom.ne.'P '.and. &
4067 atom(1:1).ne.'H' .and. &
4068 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
4070 ! write (iout,*) "sidechain ",atom
4071 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4072 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
4073 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4075 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4079 ! print *,"IONS",ions,card(1:6)
4080 else if ((ions).and.(card(1:6).eq.'HETATM')) then
4081 if (firstion.eq.0) then
4085 dc(j,ires)=sccor(j,iii)
4088 call sccenter(ires,iii,sccor)
4091 read (card(12:16),*) atom
4092 print *,"HETATOM", atom(1:2)
4093 read (card(18:20),'(a3)') res
4094 if (atom(3:3).eq.'H') cycle
4095 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
4096 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
4097 .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
4098 (atom(1:2).eq.'O ')) then
4100 print *,"I have water"
4101 if (molecule.ne.5) molecprev=molecule
4103 nres_molec(molecule)=nres_molec(molecule)+1
4104 print *,"HERE",nres_molec(molecule)
4105 if (res.ne.'WAT') res=res(2:3)//' '
4106 itype(ires,molecule)=rescode(ires,res,0,molecule)
4107 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4111 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
4112 if (ires.eq.0) return
4113 ! Calculate dummy residue coordinates inside the "chain" of a multichain
4116 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
4118 nres_molec(molecule)=nres_molec(molecule)-2
4119 print *,'I have',nres, nres_molec(:)
4121 do k=1,4 ! ions are without dummy
4122 if (nres_molec(k).eq.0) cycle
4124 ! write (iout,*) i,itype(i,1)
4125 ! if (itype(i,1).eq.ntyp1) then
4126 ! write (iout,*) "dummy",i,itype(i,1)
4128 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
4129 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
4133 if (itype(i,k).eq.ntyp1_molec(k)) then
4134 if (itype(i+1,k).eq.ntyp1_molec(k)) then
4135 if (itype(i+2,k).eq.0) then
4136 ! print *,"masz sieczke"
4138 if (itype(i+2,j).ne.0) then
4140 itype(i+1,j)=ntyp1_molec(j)
4141 nres_molec(k)=nres_molec(k)-1
4142 nres_molec(j)=nres_molec(j)+1
4148 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
4149 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
4150 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
4151 ! if (unres_pdb) then
4152 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4153 ! print *,i,'tu dochodze'
4154 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
4162 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
4166 dcj=(c(j,i-2)-c(j,i-3))/2.0
4167 if (dcj.eq.0) dcj=1.23591524223
4172 else !itype(i+1,1).eq.ntyp1
4173 ! if (unres_pdb) then
4174 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4175 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
4182 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
4183 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
4187 dcj=(c(j,i+3)-c(j,i+2))/2.0
4188 if (dcj.eq.0) dcj=1.23591524223
4193 endif !itype(i+1,1).eq.ntyp1
4194 endif !itype.eq.ntyp1
4198 ! Calculate the CM of the last side chain.
4202 dc(j,ires)=sccor(j,iii)
4205 call sccenter(ires,iii,sccor)
4211 ! print *,"molecule",molecule
4212 if ((itype(nres,1).ne.10)) then
4215 if (molecule.eq.5) molecule=molecprev
4216 itype(nres,molecule)=ntyp1_molec(molecule)
4217 nres_molec(molecule)=nres_molec(molecule)+1
4218 ! if (unres_pdb) then
4219 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4220 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
4227 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
4231 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
4232 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
4233 c(j,2*nres)=c(j,nres)
4237 ! print *,'I have',nres, nres_molec(:)
4239 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
4241 if (nres.ne.nres0) then
4242 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
4244 stop "Error nres value in WHAM input"
4247 !---------------------------------
4248 !el reallocate tables
4251 ! hfrag_alloc(j,i)=hfrag(j,i)
4254 ! bfrag_alloc(j,i)=bfrag(j,i)
4260 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
4261 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
4262 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
4263 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
4267 ! hfrag(j,i)=hfrag_alloc(j,i)
4272 ! bfrag(j,i)=bfrag_alloc(j,i)
4275 !el end reallocate tables
4276 !---------------------------------
4284 c(j,2*nres)=c(j,nres)
4287 if (itype(1,1).eq.ntyp1) then
4291 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4292 call refsys(2,3,4,e1,e2,e3,fail)
4299 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4300 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4304 dcj=(c(j,4)-c(j,3))/2.0
4310 ! First lets assign correct dummy to correct type of chain
4312 if (itype(1,1).eq.ntyp1) then
4313 if (itype(2,1).eq.0) then
4315 if (itype(2,j).ne.0) then
4317 itype(1,j)=ntyp1_molec(j)
4318 nres_molec(1)=nres_molec(1)-1
4319 nres_molec(j)=nres_molec(j)+1
4326 print *,'I have',nres, nres_molec(:)
4328 ! Copy the coordinates to reference coordinates
4334 ! Calculate internal coordinates.
4336 write (iout,'(/a)') &
4337 "Cartesian coordinates of the reference structure"
4338 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4339 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4341 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4342 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4343 (c(j,ires+nres),j=1,3)
4346 ! znamy już nres wiec mozna alokowac tablice
4347 ! Calculate internal coordinates.
4348 if(me.eq.king.or..not.out1file)then
4349 write (iout,'(a)') &
4350 "Backbone and SC coordinates as read from the PDB"
4352 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4353 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4354 (c(j,nres+ires),j=1,3)
4357 ! NOW LETS ROCK! SORTING
4358 allocate(c_temporary(3,2*nres))
4359 allocate(itype_temporary(nres,5))
4360 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4361 if (.not.allocated(istype)) write(iout,*) &
4362 "SOMETHING WRONG WITH ISTYTPE"
4363 allocate(istype_temp(nres))
4364 itype_temporary(:,:)=0
4368 if (itype(i,k).ne.0) then
4370 c_temporary(j,seqalingbegin)=c(j,i)
4371 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4373 itype_temporary(seqalingbegin,k)=itype(i,k)
4374 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4375 istype_temp(seqalingbegin)=istype(i)
4376 molnum(seqalingbegin)=k
4377 seqalingbegin=seqalingbegin+1
4383 c(j,i)=c_temporary(j,i)
4385 if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
4391 itype(i,k)=itype_temporary(i,k)
4392 istype(i)=istype_temp(i)
4395 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4396 ! I have only ions now dummy atoms in the system
4398 itype(1,5)=itype(2,5)
4404 itype(i,5)=itype(i+1,5)
4411 nres_molec(1)=nres_molec(1)-1
4413 ! if (itype(1,1).eq.ntyp1) then
4416 ! if (unres_pdb) then
4417 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4418 ! call refsys(2,3,4,e1,e2,e3,fail)
4425 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4429 ! dcj=(c(j,4)-c(j,3))/2.0
4431 ! c(j,nres+1)=c(j,1)
4437 write (iout,'(/a)') &
4438 "Cartesian coordinates of the reference structure after sorting"
4439 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4440 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4442 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4443 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4444 (c(j,ires+nres),j=1,3)
4448 print *,seqalingbegin,nres
4449 if(.not.allocated(vbld)) then
4450 allocate(vbld(2*nres))
4455 if(.not.allocated(vbld_inv)) then
4456 allocate(vbld_inv(2*nres))
4462 if(.not.allocated(theta)) then
4463 allocate(theta(nres+2))
4467 if(.not.allocated(phi)) allocate(phi(nres+2))
4468 if(.not.allocated(alph)) allocate(alph(nres+2))
4469 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4470 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4471 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4472 if(.not.allocated(costtab)) allocate(costtab(nres))
4473 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4474 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4475 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4476 if(.not.allocated(xxref)) allocate(xxref(nres))
4477 if(.not.allocated(yyref)) allocate(yyref(nres))
4478 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4479 if(.not.allocated(dc_norm)) then
4480 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4481 allocate(dc_norm(3,0:2*nres+2))
4484 write(iout,*) "before int_from_cart"
4485 call int_from_cart(.true.,.false.)
4486 call sc_loc_geom(.false.)
4487 write(iout,*) "after int_from_cart"
4491 thetaref(i)=theta(i)
4494 write(iout,*) "after thetaref"
4502 dc(j,i)=c(j,i+1)-c(j,i)
4503 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4508 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4509 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4511 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4515 ! Copy the coordinates to reference coordinates
4516 ! Splits to single chain if occurs
4520 ! cref(j,i,cou)=c(j,i)
4526 ! chomo(j,i,k)=c(j,i)
4529 ! write(iout,*) "after chomo"
4531 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4532 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4533 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4534 !-----------------------------
4538 write (iout,*) "symetr", symetr
4541 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4543 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4546 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4551 cref(j,i,cou)=c(j,i)
4552 cref(j,i+nres,cou)=c(j,i+nres)
4554 chain_rep(j,lll,kkk)=c(j,i)
4555 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4559 write (iout,*) chain_length
4560 if (chain_length.eq.0) chain_length=nres
4562 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4563 chain_rep(j,chain_length+nres,symetr) &
4564 =chain_rep(j,chain_length+nres,1)
4567 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4569 ! do kkk=1,chain_length
4570 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4574 ! makes copy of chains
4575 write (iout,*) "symetr", symetr
4580 if (symetr.gt.1) then
4587 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4593 write (iout,*) i,icha
4594 do lll=1,chain_length
4596 if (cou.le.nres) then
4598 kupa=mod(lll,chain_length)
4599 iprzes=(kkk-1)*chain_length+lll
4600 if (kupa.eq.0) kupa=chain_length
4601 write (iout,*) "kupa", kupa
4602 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4603 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4610 !-koniec robienia kopii
4613 write (iout,*) "nowa struktura", nperm
4615 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4617 cref(3,i,kkk),cref(1,nres+i,kkk),&
4618 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4620 100 format (//' alpha-carbon coordinates ',&
4621 ' centroid coordinates'/ &
4622 ' ', 6X,'X',11X,'Y',11X,'Z', &
4623 10X,'X',11X,'Y',11X,'Z')
4624 110 format (a,'(',i5,')',6f12.5)
4630 bfrag(i,j)=bfrag(i,j)-ishift
4636 hfrag(i,j)=hfrag(i,j)-ishift
4642 end subroutine readpdb
4643 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4644 !-----------------------------------------------------------------------------
4646 !-----------------------------------------------------------------------------
4647 subroutine read_control
4661 use random, only: random_init
4662 ! implicit real*8 (a-h,o-z)
4663 ! include 'DIMENSIONS'
4665 use prng, only:prng_restart
4667 logical :: OKRandom!, prng_restart
4670 ! include 'COMMON.IOUNITS'
4671 ! include 'COMMON.TIME1'
4672 ! include 'COMMON.THREAD'
4673 ! include 'COMMON.SBRIDGE'
4674 ! include 'COMMON.CONTROL'
4675 ! include 'COMMON.MCM'
4676 ! include 'COMMON.MAP'
4677 ! include 'COMMON.HEADER'
4678 ! include 'COMMON.CSA'
4679 ! include 'COMMON.CHAIN'
4680 ! include 'COMMON.MUCA'
4681 ! include 'COMMON.MD'
4682 ! include 'COMMON.FFIELD'
4683 ! include 'COMMON.INTERACT'
4684 ! include 'COMMON.SETUP'
4685 !el integer :: KDIAG,ICORFL,IXDR
4686 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4687 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4688 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4689 ! character(len=80) :: ucase
4690 character(len=640) :: controlcard
4692 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4694 usampl=.false. ! the default value of usample should be 0
4695 ! write(iout,*) "KURWA2",usampl
4699 read (INP,'(a)') titel
4700 call card_concat(controlcard,.true.)
4701 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4702 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4703 call reada(controlcard,'SEED',seed,0.0D0)
4704 call random_init(seed)
4705 ! Set up the time limit (caution! The time must be input in minutes!)
4706 read_cart=index(controlcard,'READ_CART').gt.0
4707 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4708 call readi(controlcard,'SYM',symetr,1)
4709 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4710 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4711 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4712 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4713 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4714 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4715 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4716 call reada(controlcard,'DRMS',drms,0.1D0)
4717 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4718 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4719 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4720 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4721 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4722 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4723 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4724 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4725 write (iout,'(a,f10.1)')'DRMS = ',drms
4726 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4727 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4729 call readi(controlcard,'NZ_START',nz_start,0)
4730 call readi(controlcard,'NZ_END',nz_end,0)
4731 call readi(controlcard,'IZ_SC',iz_sc,0)
4732 timlim=60.0D0*timlim
4733 safety = 60.0d0*safety
4736 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4737 !C SHIELD keyword sets if the shielding effect of side-chains is used
4738 !C 0 denots no shielding is used all peptide are equally despite the
4739 !C solvent accesible area
4740 !C 1 the newly introduced function
4741 !C 2 reseved for further possible developement
4742 call readi(controlcard,'SHIELD',shield_mode,0)
4743 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4744 write(iout,*) "shield_mode",shield_mode
4745 call readi(controlcard,'TORMODE',tor_mode,0)
4746 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4747 write(iout,*) "torsional and valence angle mode",tor_mode
4749 !C Varibles set size of box
4750 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4751 protein=index(controlcard,"PROTEIN").gt.0
4752 ions=index(controlcard,"IONS").gt.0
4753 fodson=index(controlcard,"FODSON").gt.0
4754 call readi(controlcard,'OLDION',oldion,1)
4755 nucleic=index(controlcard,"NUCLEIC").gt.0
4756 write (iout,*) "with_theta_constr ",with_theta_constr
4757 AFMlog=(index(controlcard,'AFM'))
4758 selfguide=(index(controlcard,'SELFGUIDE'))
4759 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4760 call readi(controlcard,'GENCONSTR',genconstr,0)
4761 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4762 call reada(controlcard,'BOXY',boxysize,100.0d0)
4763 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4764 call readi(controlcard,'TUBEMOD',tubemode,0)
4765 print *,"SCELE",scelemode
4766 call readi(controlcard,"SCELEMODE",scelemode,0)
4767 print *,"SCELE",scelemode
4769 ! elemode = 0 is orignal UNRES electrostatics
4770 ! elemode = 1 is "Momo" potentials in progress
4771 ! elemode = 2 is in development EVALD
4774 write (iout,*) TUBEmode,"TUBEMODE"
4775 if (TUBEmode.gt.0) then
4776 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4777 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4778 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4779 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4780 call reada(controlcard,"VNANO",velnanoconst,0.0d0)
4781 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4782 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4783 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4784 buftubebot=bordtubebot+tubebufthick
4785 buftubetop=bordtubetop-tubebufthick
4788 ! CUTOFFF ON ELECTROSTATICS
4789 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4790 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4791 write(iout,*) "R_CUT_ELE=",r_cut_ele
4792 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
4793 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
4794 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
4796 ! Lipidec parameters
4797 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4798 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4799 if (lipthick.gt.0.0d0) then
4800 bordliptop=(boxzsize+lipthick)/2.0
4801 bordlipbot=bordliptop-lipthick
4802 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4803 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4804 buflipbot=bordlipbot+lipbufthick
4805 bufliptop=bordliptop-lipbufthick
4806 if ((lipbufthick*2.0d0).gt.lipthick) &
4807 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4808 endif !lipthick.gt.0
4809 write(iout,*) "bordliptop=",bordliptop
4810 write(iout,*) "bordlipbot=",bordlipbot
4811 write(iout,*) "bufliptop=",bufliptop
4812 write(iout,*) "buflipbot=",buflipbot
4813 write (iout,*) "SHIELD MODE",shield_mode
4815 !C-------------------------
4816 minim=(index(controlcard,'MINIMIZE').gt.0)
4817 dccart=(index(controlcard,'CART').gt.0)
4818 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4819 overlapsc=.not.overlapsc
4820 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4821 searchsc=.not.searchsc
4822 sideadd=(index(controlcard,'SIDEADD').gt.0)
4823 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4824 outpdb=(index(controlcard,'PDBOUT').gt.0)
4825 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4826 pdbref=(index(controlcard,'PDBREF').gt.0)
4827 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4828 indpdb=index(controlcard,'PDBSTART')
4829 extconf=(index(controlcard,'EXTCONF').gt.0)
4830 call readi(controlcard,'IPRINT',iprint,0)
4831 call readi(controlcard,'MAXGEN',maxgen,10000)
4832 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4833 call readi(controlcard,"KDIAG",kdiag,0)
4834 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4835 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4836 write (iout,*) "RESCALE_MODE",rescale_mode
4837 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4838 if (index(controlcard,'REGULAR').gt.0.0D0) then
4839 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4843 if (index(controlcard,'CHECKGRAD').gt.0) then
4845 if (index(controlcard,'CART').gt.0) then
4847 elseif (index(controlcard,'CARINT').gt.0) then
4852 elseif (index(controlcard,'THREAD').gt.0) then
4854 call readi(controlcard,'THREAD',nthread,0)
4855 if (nthread.gt.0) then
4856 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4859 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4860 stop 'Error termination in Read_Control.'
4862 else if (index(controlcard,'MCMA').gt.0) then
4864 else if (index(controlcard,'MCEE').gt.0) then
4866 else if (index(controlcard,'MULTCONF').gt.0) then
4868 else if (index(controlcard,'MAP').gt.0) then
4870 call readi(controlcard,'MAP',nmap,0)
4871 else if (index(controlcard,'CSA').gt.0) then
4873 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4875 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4878 !fcm else if (index(controlcard,'MCMF').gt.0) then
4880 else if (index(controlcard,'SOFTREG').gt.0) then
4882 else if (index(controlcard,'CHECK_BOND').gt.0) then
4884 else if (index(controlcard,'TEST').gt.0) then
4886 else if (index(controlcard,'MD').gt.0) then
4888 else if (index(controlcard,'RE ').gt.0) then
4892 lmuca=index(controlcard,'MUCA').gt.0
4893 call readi(controlcard,'MUCADYN',mucadyn,0)
4894 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4895 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4897 write (iout,*) 'MUCADYN=',mucadyn
4898 write (iout,*) 'MUCASMOOTH=',muca_smooth
4901 iscode=index(controlcard,'ONE_LETTER')
4902 indphi=index(controlcard,'PHI')
4903 indback=index(controlcard,'BACK')
4904 iranconf=index(controlcard,'RAND_CONF')
4905 i2ndstr=index(controlcard,'USE_SEC_PRED')
4906 gradout=index(controlcard,'GRADOUT').gt.0
4907 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4908 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4909 if (me.eq.king .or. .not.out1file ) &
4910 write (iout,*) "DISTCHAINMAX",distchainmax
4912 if(me.eq.king.or..not.out1file) &
4913 write (iout,'(2a)') diagmeth(kdiag),&
4914 ' routine used to diagonalize matrices.'
4915 if (shield_mode.gt.0) then
4916 pi=4.0D0*datan(1.0D0)
4917 !C VSolvSphere the volume of solving sphere
4919 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4920 !C there will be no distinction between proline peptide group and normal peptide
4921 !C group in case of shielding parameters
4922 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4923 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4924 write (iout,*) VSolvSphere,VSolvSphere_div
4925 !C long axis of side chain
4927 ! long_r_sidechain(i)=vbldsc0(1,i)
4928 ! short_r_sidechain(i)=sigma0(i)
4929 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4936 end subroutine read_control
4937 !-----------------------------------------------------------------------------
4938 subroutine read_REMDpar
4940 ! Read REMD settings
4947 use control_data, only:out1file
4948 ! implicit real*8 (a-h,o-z)
4949 ! include 'DIMENSIONS'
4950 ! include 'COMMON.IOUNITS'
4951 ! include 'COMMON.TIME1'
4952 ! include 'COMMON.MD'
4955 !el include 'COMMON.LANGEVIN'
4957 !el include 'COMMON.LANGEVIN.lang0'
4959 ! include 'COMMON.INTERACT'
4960 ! include 'COMMON.NAMES'
4961 ! include 'COMMON.GEO'
4962 ! include 'COMMON.REMD'
4963 ! include 'COMMON.CONTROL'
4964 ! include 'COMMON.SETUP'
4965 ! character(len=80) :: ucase
4966 character(len=320) :: controlcard
4967 character(len=3200) :: controlcard1
4968 integer :: iremd_m_total
4971 ! real(kind=8) :: var,ene
4973 if(me.eq.king.or..not.out1file) &
4974 write (iout,*) "REMD setup"
4976 call card_concat(controlcard,.true.)
4977 call readi(controlcard,"NREP",nrep,3)
4978 call readi(controlcard,"NSTEX",nstex,1000)
4979 call reada(controlcard,"RETMIN",retmin,10.0d0)
4980 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4981 mremdsync=(index(controlcard,'SYNC').gt.0)
4982 call readi(controlcard,"NSYN",i_sync_step,100)
4983 restart1file=(index(controlcard,'REST1FILE').gt.0)
4984 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4985 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4986 if(max_cache_traj_use.gt.max_cache_traj) &
4987 max_cache_traj_use=max_cache_traj
4988 if(me.eq.king.or..not.out1file) then
4989 !d if (traj1file) then
4990 !rc caching is in testing - NTWX is not ignored
4991 !d write (iout,*) "NTWX value is ignored"
4992 !d write (iout,*) " trajectory is stored to one file by master"
4993 !d write (iout,*) " before exchange at NSTEX intervals"
4995 write (iout,*) "NREP= ",nrep
4996 write (iout,*) "NSTEX= ",nstex
4997 write (iout,*) "SYNC= ",mremdsync
4998 write (iout,*) "NSYN= ",i_sync_step
4999 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
5002 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
5003 if (index(controlcard,'TLIST').gt.0) then
5005 call card_concat(controlcard1,.true.)
5006 read(controlcard1,*) (remd_t(i),i=1,nrep)
5007 if(me.eq.king.or..not.out1file) &
5008 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
5011 if (index(controlcard,'MLIST').gt.0) then
5013 call card_concat(controlcard1,.true.)
5014 read(controlcard1,*) (remd_m(i),i=1,nrep)
5015 if(me.eq.king.or..not.out1file) then
5016 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
5019 iremd_m_total=iremd_m_total+remd_m(i)
5021 write (iout,*) 'Total number of replicas ',iremd_m_total
5024 if(me.eq.king.or..not.out1file) &
5025 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
5027 end subroutine read_REMDpar
5028 !-----------------------------------------------------------------------------
5029 subroutine read_MDpar
5033 use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
5035 use geometry_data, only: pi
5037 ! implicit real*8 (a-h,o-z)
5038 ! include 'DIMENSIONS'
5039 ! include 'COMMON.IOUNITS'
5040 ! include 'COMMON.TIME1'
5041 ! include 'COMMON.MD'
5044 !el include 'COMMON.LANGEVIN'
5046 !el include 'COMMON.LANGEVIN.lang0'
5048 ! include 'COMMON.INTERACT'
5049 ! include 'COMMON.NAMES'
5050 ! include 'COMMON.GEO'
5051 ! include 'COMMON.SETUP'
5052 ! include 'COMMON.CONTROL'
5053 ! include 'COMMON.SPLITELE'
5054 ! character(len=80) :: ucase
5055 character(len=320) :: controlcard
5060 call card_concat(controlcard,.true.)
5061 call readi(controlcard,"NSTEP",n_timestep,1000000)
5062 call readi(controlcard,"NTWE",ntwe,100)
5063 call readi(controlcard,"NTWX",ntwx,1000)
5064 call readi(controlcard,"NFOD",nfodstep,100)
5065 call reada(controlcard,"DT",d_time,1.0d-1)
5066 call reada(controlcard,"DVMAX",dvmax,2.0d1)
5067 call reada(controlcard,"DAMAX",damax,1.0d1)
5068 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
5069 call readi(controlcard,"LANG",lang,0)
5070 RESPA = index(controlcard,"RESPA") .gt. 0
5071 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
5072 ntime_split0=ntime_split
5073 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
5074 ntime_split0=ntime_split
5075 call reada(controlcard,"R_CUT",r_cut,2.0d0)
5076 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
5077 rest = index(controlcard,"REST").gt.0
5078 tbf = index(controlcard,"TBF").gt.0
5079 usampl = index(controlcard,"USAMPL").gt.0
5080 ! write(iout,*) "KURWA",usampl
5081 mdpdb = index(controlcard,"MDPDB").gt.0
5082 call reada(controlcard,"T_BATH",t_bath,300.0d0)
5083 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
5084 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
5085 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
5086 if (count_reset_moment.eq.0) count_reset_moment=1000000000
5087 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
5088 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
5089 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
5090 if (count_reset_vel.eq.0) count_reset_vel=1000000000
5091 large = index(controlcard,"LARGE").gt.0
5092 print_compon = index(controlcard,"PRINT_COMPON").gt.0
5093 rattle = index(controlcard,"RATTLE").gt.0
5094 preminim=(index(controlcard,'PREMINIM').gt.0)
5095 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
5096 write (iout,*) "PREMINIM ",preminim
5097 dccart=(index(controlcard,'CART').gt.0)
5098 if (preminim) call read_minim
5099 ! if performing umbrella sampling, fragments constrained are read from the fragment file
5105 if(me.eq.king.or..not.out1file) then
5107 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
5109 write (iout,'(a)') "The units are:"
5110 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
5111 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
5112 " acceleration: angstrom/(48.9 fs)**2"
5113 write (iout,'(a)') "energy: kcal/mol, temperature: K"
5115 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
5116 write (iout,'(a60,f10.5,a)') &
5117 "Initial time step of numerical integration:",d_time,&
5119 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
5121 write (iout,'(2a,i4,a)') &
5122 "A-MTS algorithm used; initial time step for fast-varying",&
5123 " short-range forces split into",ntime_split," steps."
5124 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
5125 r_cut," lambda",rlamb
5127 write (iout,'(2a,f10.5)') &
5128 "Maximum acceleration threshold to reduce the time step",&
5129 "/increase split number:",damax
5130 write (iout,'(2a,f10.5)') &
5131 "Maximum predicted energy drift to reduce the timestep",&
5132 "/increase split number:",edriftmax
5133 write (iout,'(a60,f10.5)') &
5134 "Maximum velocity threshold to reduce velocities:",dvmax
5135 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
5136 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
5137 if (rattle) write (iout,'(a60)') &
5138 "Rattle algorithm used to constrain the virtual bonds"
5142 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
5143 call reada(controlcard,"RWAT",rwat,1.4d0)
5144 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
5145 surfarea=index(controlcard,"SURFAREA").gt.0
5146 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
5147 if(me.eq.king.or..not.out1file)then
5148 write (iout,'(/a,$)') "Langevin dynamics calculation"
5150 write (iout,'(a/)') &
5151 " with direct integration of Langevin equations"
5152 else if (lang.eq.2) then
5153 write (iout,'(a/)') " with TINKER stochasic MD integrator"
5154 else if (lang.eq.3) then
5155 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
5156 else if (lang.eq.4) then
5157 write (iout,'(a/)') " in overdamped mode"
5159 write (iout,'(//a,i5)') &
5160 "=========== ERROR: Unknown Langevin dynamics mode:",lang
5163 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5164 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
5165 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
5166 write (iout,'(a60,f10.5)') &
5167 "Scaling factor of the friction forces:",scal_fric
5168 if (surfarea) write (iout,'(2a,i10,a)') &
5169 "Friction coefficients will be scaled by solvent-accessible",&
5170 " surface area every",reset_fricmat," steps."
5172 ! Calculate friction coefficients and bounds of stochastic forces
5173 eta=6*pi*cPoise*etawat
5174 if(me.eq.king.or..not.out1file) &
5175 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
5178 do j=1,5 !types of molecules
5179 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
5180 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
5182 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
5183 do j=1,5 !types of molecules
5185 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
5186 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
5190 if(me.eq.king.or..not.out1file)then
5192 write (iout,'(/2a/)') &
5193 "Radii of site types and friction coefficients and std's of",&
5194 " stochastic forces of fully exposed sites"
5195 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
5198 write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
5199 gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
5204 if(me.eq.king.or..not.out1file)then
5205 write (iout,'(a)') "Berendsen bath calculation"
5206 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5207 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
5209 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
5210 count_reset_moment," steps"
5212 write (iout,'(a,i10,a)') &
5213 "Velocities will be reset at random every",count_reset_vel,&
5217 if(me.eq.king.or..not.out1file) &
5218 write (iout,'(a31)') "Microcanonical mode calculation"
5220 if(me.eq.king.or..not.out1file)then
5221 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
5223 write(iout,*) "MD running with constraints."
5224 write(iout,*) "Equilibration time ", eq_time, " mtus."
5225 write(iout,*) "Constraining ", nfrag," fragments."
5226 write(iout,*) "Length of each fragment, weight and q0:"
5228 write (iout,*) "Set of restraints #",iset
5230 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
5231 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
5233 write(iout,*) "constraints between ", npair, "fragments."
5234 write(iout,*) "constraint pairs, weights and q0:"
5236 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
5237 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
5239 write(iout,*) "angle constraints within ", nfrag_back,&
5240 "backbone fragments."
5241 write(iout,*) "fragment, weights:"
5243 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
5244 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
5245 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
5248 iset=mod(kolor,nset)+1
5251 if(me.eq.king.or..not.out1file) &
5252 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
5254 end subroutine read_MDpar
5255 !-----------------------------------------------------------------------------
5259 ! implicit real*8 (a-h,o-z)
5260 ! include 'DIMENSIONS'
5261 ! include 'COMMON.MAP'
5262 ! include 'COMMON.IOUNITS'
5263 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
5264 character(len=80) :: mapcard !,ucase
5267 ! real(kind=8) :: var,ene
5270 read (inp,'(a)') mapcard
5271 mapcard=ucase(mapcard)
5272 if (index(mapcard,'PHI').gt.0) then
5274 else if (index(mapcard,'THE').gt.0) then
5276 else if (index(mapcard,'ALP').gt.0) then
5278 else if (index(mapcard,'OME').gt.0) then
5281 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
5282 stop 'Error - illegal variable spec in MAP card.'
5284 call readi (mapcard,'RES1',res1(imap),0)
5285 call readi (mapcard,'RES2',res2(imap),0)
5286 if (res1(imap).eq.0) then
5287 res1(imap)=res2(imap)
5288 else if (res2(imap).eq.0) then
5289 res2(imap)=res1(imap)
5291 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
5292 write (iout,'(a)') &
5293 'Error - illegal definition of variable group in MAP.'
5294 stop 'Error - illegal definition of variable group in MAP.'
5296 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
5297 call reada(mapcard,'TO',ang_to(imap),0.0D0)
5298 call readi(mapcard,'NSTEP',nstep(imap),0)
5299 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
5300 write (iout,'(a)') &
5301 'Illegal boundary and/or step size specification in MAP.'
5302 stop 'Illegal boundary and/or step size specification in MAP.'
5306 end subroutine map_read
5307 !-----------------------------------------------------------------------------
5310 use control_data, only: vdisulf
5312 ! implicit real*8 (a-h,o-z)
5313 ! include 'DIMENSIONS'
5314 ! include 'COMMON.IOUNITS'
5315 ! include 'COMMON.GEO'
5316 ! include 'COMMON.CSA'
5317 ! include 'COMMON.BANK'
5318 ! include 'COMMON.CONTROL'
5319 ! character(len=80) :: ucase
5320 character(len=620) :: mcmcard
5322 ! integer :: ntf,ik,iw_pdb
5323 ! real(kind=8) :: var,ene
5325 call card_concat(mcmcard,.true.)
5327 call readi(mcmcard,'NCONF',nconf,50)
5328 call readi(mcmcard,'NADD',nadd,0)
5329 call readi(mcmcard,'JSTART',jstart,1)
5330 call readi(mcmcard,'JEND',jend,1)
5331 call readi(mcmcard,'NSTMAX',nstmax,500000)
5332 call readi(mcmcard,'N0',n0,1)
5333 call readi(mcmcard,'N1',n1,6)
5334 call readi(mcmcard,'N2',n2,4)
5335 call readi(mcmcard,'N3',n3,0)
5336 call readi(mcmcard,'N4',n4,0)
5337 call readi(mcmcard,'N5',n5,0)
5338 call readi(mcmcard,'N6',n6,10)
5339 call readi(mcmcard,'N7',n7,0)
5340 call readi(mcmcard,'N8',n8,0)
5341 call readi(mcmcard,'N9',n9,0)
5342 call readi(mcmcard,'N14',n14,0)
5343 call readi(mcmcard,'N15',n15,0)
5344 call readi(mcmcard,'N16',n16,0)
5345 call readi(mcmcard,'N17',n17,0)
5346 call readi(mcmcard,'N18',n18,0)
5348 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5350 call readi(mcmcard,'NDIFF',ndiff,2)
5351 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5352 call readi(mcmcard,'IS1',is1,1)
5353 call readi(mcmcard,'IS2',is2,8)
5354 call readi(mcmcard,'NRAN0',nran0,4)
5355 call readi(mcmcard,'NRAN1',nran1,2)
5356 call readi(mcmcard,'IRR',irr,1)
5357 call readi(mcmcard,'NSEED',nseed,20)
5358 call readi(mcmcard,'NTOTAL',ntotal,10000)
5359 call reada(mcmcard,'CUT1',cut1,2.0d0)
5360 call reada(mcmcard,'CUT2',cut2,5.0d0)
5361 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5362 call readi(mcmcard,'ICMAX',icmax,3)
5363 call readi(mcmcard,'IRESTART',irestart,0)
5364 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5367 call reada(mcmcard,'DELE',dele,20.0d0)
5368 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5369 call readi(mcmcard,'IREF',iref,0)
5370 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5371 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5372 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5373 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5374 write (iout,*) "NCONF_IN",nconf_in
5376 end subroutine csaread
5377 !-----------------------------------------------------------------------------
5381 use control_data, only: MaxMoveType
5384 ! implicit real*8 (a-h,o-z)
5385 ! include 'DIMENSIONS'
5386 ! include 'COMMON.MCM'
5387 ! include 'COMMON.MCE'
5388 ! include 'COMMON.IOUNITS'
5389 ! character(len=80) :: ucase
5390 character(len=320) :: mcmcard
5393 ! real(kind=8) :: var,ene
5395 call card_concat(mcmcard,.true.)
5396 call readi(mcmcard,'MAXACC',maxacc,100)
5397 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5398 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5399 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5400 call readi(mcmcard,'MAXREPM',maxrepm,200)
5401 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5402 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5403 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5404 call reada(mcmcard,'E_UP',e_up,5.0D0)
5405 call reada(mcmcard,'DELTE',delte,0.1D0)
5406 call readi(mcmcard,'NSWEEP',nsweep,5)
5407 call readi(mcmcard,'NSTEPH',nsteph,0)
5408 call readi(mcmcard,'NSTEPC',nstepc,0)
5409 call reada(mcmcard,'TMIN',tmin,298.0D0)
5410 call reada(mcmcard,'TMAX',tmax,298.0D0)
5411 call readi(mcmcard,'NWINDOW',nwindow,0)
5412 call readi(mcmcard,'PRINT_MC',print_mc,0)
5413 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5414 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5415 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5416 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5417 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5418 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5419 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5420 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5421 if (nwindow.gt.0) then
5422 allocate(winstart(nwindow)) !!el (maxres)
5423 allocate(winend(nwindow)) !!el
5424 allocate(winlen(nwindow)) !!el
5425 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5427 winlen(i)=winend(i)-winstart(i)+1
5430 if (tmax.lt.tmin) tmax=tmin
5431 if (tmax.eq.tmin) then
5435 if (nstepc.gt.0 .and. nsteph.gt.0) then
5436 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5437 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5439 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5440 ! Probabilities of different move types
5441 sumpro_type(0)=0.0D0
5442 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5443 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5444 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5445 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5446 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5447 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5448 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5450 print *,'i',i,' sumprotype',sumpro_type(i)
5451 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5452 print *,'i',i,' sumprotype',sumpro_type(i)
5455 end subroutine mcmread
5456 !-----------------------------------------------------------------------------
5457 subroutine read_minim
5460 ! implicit real*8 (a-h,o-z)
5461 ! include 'DIMENSIONS'
5462 ! include 'COMMON.MINIM'
5463 ! include 'COMMON.IOUNITS'
5464 ! character(len=80) :: ucase
5465 character(len=320) :: minimcard
5467 ! integer :: ntf,ik,iw_pdb
5468 ! real(kind=8) :: var,ene
5470 call card_concat(minimcard,.true.)
5471 call readi(minimcard,'MAXMIN',maxmin,2000)
5472 call readi(minimcard,'MAXFUN',maxfun,5000)
5473 call readi(minimcard,'MINMIN',minmin,maxmin)
5474 call readi(minimcard,'MINFUN',minfun,maxmin)
5475 call reada(minimcard,'TOLF',tolf,1.0D-2)
5476 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5477 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5478 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5479 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5480 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5481 'Options in energy minimization:'
5482 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5483 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5484 'MinMin:',MinMin,' MinFun:',MinFun,&
5485 ' TolF:',TolF,' RTolF:',RTolF
5487 end subroutine read_minim
5488 !-----------------------------------------------------------------------------
5489 subroutine openunits
5491 use MD_data, only: usampl
5494 use control_data, only:out1file
5495 use control, only: getenv_loc
5496 ! implicit real*8 (a-h,o-z)
5497 ! include 'DIMENSIONS'
5500 character(len=16) :: form,nodename
5501 integer :: nodelen,ierror,npos
5503 ! include 'COMMON.SETUP'
5504 ! include 'COMMON.IOUNITS'
5505 ! include 'COMMON.MD'
5506 ! include 'COMMON.CONTROL'
5507 integer :: lenpre,lenpot,lentmp !,ilen
5509 character(len=3) :: out1file_text !,ucase
5510 character(len=3) :: ll
5513 ! integer :: ntf,ik,iw_pdb
5514 ! real(kind=8) :: var,ene
5516 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5517 call getenv_loc("PREFIX",prefix)
5519 call getenv_loc("POT",pot)
5520 call getenv_loc("DIRTMP",tmpdir)
5521 call getenv_loc("CURDIR",curdir)
5522 call getenv_loc("OUT1FILE",out1file_text)
5523 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5524 out1file_text=ucase(out1file_text)
5525 if (out1file_text(1:1).eq."Y") then
5528 out1file=fg_rank.gt.0
5533 if (lentmp.gt.0) then
5534 write (*,'(80(1h!))')
5535 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5536 write (*,'(80(1h!))')
5537 write (*,*)"All output files will be on node /tmp directory."
5539 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5540 if (me.eq.king) then
5541 write (*,*) "The master node is ",nodename
5542 else if (fg_rank.eq.0) then
5543 write (*,*) "I am the CG slave node ",nodename
5545 write (*,*) "I am the FG slave node ",nodename
5548 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5549 lenpre = lentmp+lenpre+1
5551 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5552 ! Get the names and open the input files
5553 #if defined(WINIFL) || defined(WINPGI)
5554 open(1,file=pref_orig(:ilen(pref_orig))// &
5555 '.inp',status='old',readonly,shared)
5556 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5557 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5558 ! Get parameter filenames and open the parameter files.
5559 call getenv_loc('BONDPAR',bondname)
5560 open (ibond,file=bondname,status='old',readonly,shared)
5561 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5562 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5563 call getenv_loc('THETPAR',thetname)
5564 open (ithep,file=thetname,status='old',readonly,shared)
5565 call getenv_loc('ROTPAR',rotname)
5566 open (irotam,file=rotname,status='old',readonly,shared)
5567 call getenv_loc('TORPAR',torname)
5568 open (itorp,file=torname,status='old',readonly,shared)
5569 call getenv_loc('TORDPAR',tordname)
5570 open (itordp,file=tordname,status='old',readonly,shared)
5571 call getenv_loc('FOURIER',fouriername)
5572 open (ifourier,file=fouriername,status='old',readonly,shared)
5573 call getenv_loc('ELEPAR',elename)
5574 open (ielep,file=elename,status='old',readonly,shared)
5575 call getenv_loc('SIDEPAR',sidename)
5576 open (isidep,file=sidename,status='old',readonly,shared)
5578 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5579 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5580 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5581 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5582 call getenv_loc('TORPAR_NUCL',torname_nucl)
5583 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5584 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5585 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5586 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5587 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5590 #elif (defined CRAY) || (defined AIX)
5591 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5593 ! print *,"Processor",myrank," opened file 1"
5594 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5595 ! print *,"Processor",myrank," opened file 9"
5596 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5597 ! Get parameter filenames and open the parameter files.
5598 call getenv_loc('BONDPAR',bondname)
5599 open (ibond,file=bondname,status='old',action='read')
5600 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5601 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5603 ! print *,"Processor",myrank," opened file IBOND"
5604 call getenv_loc('THETPAR',thetname)
5605 open (ithep,file=thetname,status='old',action='read')
5606 ! print *,"Processor",myrank," opened file ITHEP"
5607 call getenv_loc('ROTPAR',rotname)
5608 open (irotam,file=rotname,status='old',action='read')
5609 ! print *,"Processor",myrank," opened file IROTAM"
5610 call getenv_loc('TORPAR',torname)
5611 open (itorp,file=torname,status='old',action='read')
5612 ! print *,"Processor",myrank," opened file ITORP"
5613 call getenv_loc('TORDPAR',tordname)
5614 open (itordp,file=tordname,status='old',action='read')
5615 ! print *,"Processor",myrank," opened file ITORDP"
5616 call getenv_loc('SCCORPAR',sccorname)
5617 open (isccor,file=sccorname,status='old',action='read')
5618 ! print *,"Processor",myrank," opened file ISCCOR"
5619 call getenv_loc('FOURIER',fouriername)
5620 open (ifourier,file=fouriername,status='old',action='read')
5621 ! print *,"Processor",myrank," opened file IFOURIER"
5622 call getenv_loc('ELEPAR',elename)
5623 open (ielep,file=elename,status='old',action='read')
5624 ! print *,"Processor",myrank," opened file IELEP"
5625 call getenv_loc('SIDEPAR',sidename)
5626 open (isidep,file=sidename,status='old',action='read')
5628 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5629 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5630 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5631 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5632 call getenv_loc('TORPAR_NUCL',torname_nucl)
5633 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5634 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5635 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5636 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5637 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5639 call getenv_loc('LIPTRANPAR',liptranname)
5640 open (iliptranpar,file=liptranname,status='old',action='read')
5641 call getenv_loc('TUBEPAR',tubename)
5642 open (itube,file=tubename,status='old',action='read')
5643 call getenv_loc('IONPAR',ionname)
5644 open (iion,file=ionname,status='old',action='read')
5645 call getenv_loc('IONPAR_TRAN',iontranname)
5646 open (iiontran,file=iontranname,status='old',action='read')
5647 ! print *,"Processor",myrank," opened file ISIDEP"
5648 ! print *,"Processor",myrank," opened parameter files"
5650 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5651 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5652 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5653 ! Get parameter filenames and open the parameter files.
5654 call getenv_loc('BONDPAR',bondname)
5655 open (ibond,file=bondname,status='old')
5656 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5657 open (ibond_nucl,file=bondname_nucl,status='old')
5659 call getenv_loc('THETPAR',thetname)
5660 open (ithep,file=thetname,status='old')
5661 call getenv_loc('ROTPAR',rotname)
5662 open (irotam,file=rotname,status='old')
5663 call getenv_loc('TORPAR',torname)
5664 open (itorp,file=torname,status='old')
5665 call getenv_loc('TORDPAR',tordname)
5666 open (itordp,file=tordname,status='old')
5667 call getenv_loc('SCCORPAR',sccorname)
5668 open (isccor,file=sccorname,status='old')
5669 call getenv_loc('FOURIER',fouriername)
5670 open (ifourier,file=fouriername,status='old')
5671 call getenv_loc('ELEPAR',elename)
5672 open (ielep,file=elename,status='old')
5673 call getenv_loc('SIDEPAR',sidename)
5674 open (isidep,file=sidename,status='old')
5676 open (ithep_nucl,file=thetname_nucl,status='old')
5677 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5678 open (irotam_nucl,file=rotname_nucl,status='old')
5679 call getenv_loc('TORPAR_NUCL',torname_nucl)
5680 open (itorp_nucl,file=torname_nucl,status='old')
5681 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5682 open (itordp_nucl,file=tordname_nucl,status='old')
5683 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5684 open (isidep_nucl,file=sidename_nucl,status='old')
5686 call getenv_loc('LIPTRANPAR',liptranname)
5687 open (iliptranpar,file=liptranname,status='old')
5688 call getenv_loc('TUBEPAR',tubename)
5689 open (itube,file=tubename,status='old')
5690 call getenv_loc('IONPAR',ionname)
5691 open (iion,file=ionname,status='old')
5692 call getenv_loc('IONPAR_NUCL',ionnuclname)
5693 open (iionnucl,file=ionnuclname,status='old')
5694 call getenv_loc('IONPAR_TRAN',iontranname)
5695 open (iiontran,file=iontranname,status='old')
5696 call getenv_loc('WATWAT',iwaterwatername)
5697 open (iwaterwater,file=iwaterwatername,status='old')
5698 call getenv_loc('WATPROT',iwaterscname)
5699 open (iwatersc,file=iwaterscname,status='old')
5702 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5704 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5705 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5706 ! Get parameter filenames and open the parameter files.
5707 call getenv_loc('BONDPAR',bondname)
5708 open (ibond,file=bondname,status='old',action='read')
5709 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5710 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5711 call getenv_loc('THETPAR',thetname)
5712 open (ithep,file=thetname,status='old',action='read')
5713 call getenv_loc('ROTPAR',rotname)
5714 open (irotam,file=rotname,status='old',action='read')
5715 call getenv_loc('TORPAR',torname)
5716 open (itorp,file=torname,status='old',action='read')
5717 call getenv_loc('TORDPAR',tordname)
5718 open (itordp,file=tordname,status='old',action='read')
5719 call getenv_loc('SCCORPAR',sccorname)
5720 open (isccor,file=sccorname,status='old',action='read')
5722 call getenv_loc('THETPARPDB',thetname_pdb)
5723 print *,"thetname_pdb ",thetname_pdb
5724 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5725 print *,ithep_pdb," opened"
5727 call getenv_loc('FOURIER',fouriername)
5728 open (ifourier,file=fouriername,status='old',readonly)
5729 call getenv_loc('ELEPAR',elename)
5730 open (ielep,file=elename,status='old',readonly)
5731 call getenv_loc('SIDEPAR',sidename)
5732 open (isidep,file=sidename,status='old',readonly)
5734 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5735 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5736 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5737 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5738 call getenv_loc('TORPAR_NUCL',torname_nucl)
5739 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5740 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5741 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5742 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5743 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5744 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5745 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5746 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5747 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5748 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5749 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5750 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5751 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5754 call getenv_loc('LIPTRANPAR',liptranname)
5755 open (iliptranpar,file=liptranname,status='old',action='read')
5756 call getenv_loc('LIPBOND',lipbondname)
5757 open (ilipbond,file=lipbondname,status='old',action='read')
5758 call getenv_loc('LIPNONBOND',lipnonbondname)
5759 open (ilipnonbond,file=lipnonbondname,status='old',action='read')
5760 call getenv_loc('TUBEPAR',tubename)
5761 open (itube,file=tubename,status='old',action='read')
5762 call getenv_loc('IONPAR',ionname)
5763 open (iion,file=ionname,status='old',action='read')
5764 call getenv_loc('IONPAR_NUCL',ionnuclname)
5765 open (iionnucl,file=ionnuclname,status='old',action='read')
5766 call getenv_loc('IONPAR_TRAN',iontranname)
5767 open (iiontran,file=iontranname,status='old',action='read')
5768 call getenv_loc('WATWAT',iwaterwatername)
5769 open (iwaterwater,file=iwaterwatername,status='old',action='read')
5770 call getenv_loc('WATPROT',iwaterscname)
5771 open (iwatersc,file=iwaterscname,status='old',action='read')
5774 call getenv_loc('ROTPARPDB',rotname_pdb)
5775 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5778 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5779 #if defined(WINIFL) || defined(WINPGI)
5780 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5781 #elif (defined CRAY) || (defined AIX)
5782 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5784 open (iscpp_nucl,file=scpname_nucl,status='old')
5786 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5791 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5792 ! Use -DOLDSCP to use hard-coded constants instead.
5794 call getenv_loc('SCPPAR',scpname)
5795 #if defined(WINIFL) || defined(WINPGI)
5796 open (iscpp,file=scpname,status='old',readonly,shared)
5797 #elif (defined CRAY) || (defined AIX)
5798 open (iscpp,file=scpname,status='old',action='read')
5800 open (iscpp,file=scpname,status='old')
5802 open (iscpp,file=scpname,status='old',action='read')
5805 call getenv_loc('PATTERN',patname)
5806 #if defined(WINIFL) || defined(WINPGI)
5807 open (icbase,file=patname,status='old',readonly,shared)
5808 #elif (defined CRAY) || (defined AIX)
5809 open (icbase,file=patname,status='old',action='read')
5811 open (icbase,file=patname,status='old')
5813 open (icbase,file=patname,status='old',action='read')
5816 ! Open output file only for CG processes
5817 ! print *,"Processor",myrank," fg_rank",fg_rank
5818 if (fg_rank.eq.0) then
5820 if (nodes.eq.1) then
5823 npos = dlog10(dfloat(nodes-1))+1
5825 if (npos.lt.3) npos=3
5826 write (liczba,'(i1)') npos
5827 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5829 write (liczba,form) me
5830 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5831 liczba(:ilen(liczba))
5832 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5834 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5836 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5837 liczba(:ilen(liczba))//'.mol2'
5838 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5839 liczba(:ilen(liczba))//'.stat'
5841 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5842 //liczba(:ilen(liczba))//'.stat')
5843 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5846 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5847 liczba(:ilen(liczba))//'.const'
5852 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5853 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5854 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5855 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5856 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5858 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5860 rest2name=prefix(:ilen(prefix))//'.rst'
5862 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5865 #if defined(AIX) || defined(PGI)
5866 if (me.eq.king .or. .not. out1file) &
5867 open(iout,file=outname,status='unknown')
5869 if (fg_rank.gt.0) then
5870 write (liczba,'(i3.3)') myrank/nfgtasks
5871 write (ll,'(bz,i3.3)') fg_rank
5872 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5877 open(igeom,file=intname,status='unknown',position='append')
5878 open(ipdb,file=pdbname,status='unknown')
5879 open(imol2,file=mol2name,status='unknown')
5880 open(istat,file=statname,status='unknown',position='append')
5882 !1out open(iout,file=outname,status='unknown')
5885 if (me.eq.king .or. .not.out1file) &
5886 open(iout,file=outname,status='unknown')
5888 if (fg_rank.gt.0) then
5889 write (liczba,'(i3.3)') myrank/nfgtasks
5890 write (ll,'(bz,i3.3)') fg_rank
5891 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5896 open(igeom,file=intname,status='unknown',access='append')
5897 open(ipdb,file=pdbname,status='unknown')
5898 open(imol2,file=mol2name,status='unknown')
5899 open(istat,file=statname,status='unknown',access='append')
5901 !1out open(iout,file=outname,status='unknown')
5904 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5905 csa_seed=prefix(:lenpre)//'.CSA.seed'
5906 csa_history=prefix(:lenpre)//'.CSA.history'
5907 csa_bank=prefix(:lenpre)//'.CSA.bank'
5908 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5909 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5910 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5911 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5912 csa_int=prefix(:lenpre)//'.int'
5913 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5914 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5915 csa_in=prefix(:lenpre)//'.CSA.in'
5916 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5919 write (iout,'(80(1h-))')
5920 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5921 write (iout,'(80(1h-))')
5922 write (iout,*) "Input file : ",&
5923 pref_orig(:ilen(pref_orig))//'.inp'
5924 write (iout,*) "Output file : ",&
5925 outname(:ilen(outname))
5927 write (iout,*) "Sidechain potential file : ",&
5928 sidename(:ilen(sidename))
5930 write (iout,*) "SCp potential file : ",&
5931 scpname(:ilen(scpname))
5933 write (iout,*) "Electrostatic potential file : ",&
5934 elename(:ilen(elename))
5935 write (iout,*) "Cumulant coefficient file : ",&
5936 fouriername(:ilen(fouriername))
5937 write (iout,*) "Torsional parameter file : ",&
5938 torname(:ilen(torname))
5939 write (iout,*) "Double torsional parameter file : ",&
5940 tordname(:ilen(tordname))
5941 write (iout,*) "SCCOR parameter file : ",&
5942 sccorname(:ilen(sccorname))
5943 write (iout,*) "Bond & inertia constant file : ",&
5944 bondname(:ilen(bondname))
5945 write (iout,*) "Bending parameter file : ",&
5946 thetname(:ilen(thetname))
5947 write (iout,*) "Rotamer parameter file : ",&
5948 rotname(:ilen(rotname))
5951 write (iout,*) "Thetpdb parameter file : ",&
5952 thetname_pdb(:ilen(thetname_pdb))
5955 write (iout,*) "Threading database : ",&
5956 patname(:ilen(patname))
5958 write (iout,*)" DIRTMP : ",&
5960 write (iout,'(80(1h-))')
5963 end subroutine openunits
5964 !-----------------------------------------------------------------------------
5967 use geometry_data, only: nres,dc
5969 ! implicit real*8 (a-h,o-z)
5970 ! include 'DIMENSIONS'
5971 ! include 'COMMON.CHAIN'
5972 ! include 'COMMON.IOUNITS'
5973 ! include 'COMMON.MD'
5976 ! real(kind=8) :: var,ene
5978 open(irest2,file=rest2name,status='unknown')
5979 read(irest2,*) totT,EK,potE,totE,t_bath
5982 ! AL 4/17/17: Now reading d_t(0,:) too
5984 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5987 ! AL 4/17/17: Now reading d_c(0,:) too
5989 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5992 read (irest2,*) iset
5996 end subroutine readrst
5997 !-----------------------------------------------------------------------------
5998 subroutine read_fragments
6002 use control_data, only:out1file
6005 ! implicit real*8 (a-h,o-z)
6006 ! include 'DIMENSIONS'
6010 ! include 'COMMON.SETUP'
6011 ! include 'COMMON.CHAIN'
6012 ! include 'COMMON.IOUNITS'
6013 ! include 'COMMON.MD'
6014 ! include 'COMMON.CONTROL'
6017 ! real(kind=8) :: var,ene
6019 read(inp,*) nset,nfrag,npair,nfrag_back
6021 !el from module energy
6022 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
6023 if(.not.allocated(wfrag_back)) then
6024 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6025 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6027 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
6028 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
6030 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
6031 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
6034 if(me.eq.king.or..not.out1file) &
6035 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
6036 " nfrag_back",nfrag_back
6038 read(inp,*) mset(iset)
6040 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
6042 if(me.eq.king.or..not.out1file) &
6043 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
6044 ifrag(2,i,iset), qinfrag(i,iset)
6047 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
6049 if(me.eq.king.or..not.out1file) &
6050 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
6051 ipair(2,i,iset), qinpair(i,iset)
6054 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6055 wfrag_back(3,i,iset),&
6056 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6057 if(me.eq.king.or..not.out1file) &
6058 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6059 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6063 end subroutine read_fragments
6064 !-----------------------------------------------------------------------------
6066 !-----------------------------------------------------------------------------
6070 ! implicit real*8 (a-h,o-z)
6071 ! include 'DIMENSIONS'
6072 ! include 'COMMON.CSA'
6073 ! include 'COMMON.BANK'
6074 ! include 'COMMON.IOUNITS'
6076 ! integer :: ntf,ik,iw_pdb
6077 ! real(kind=8) :: var,ene
6079 open(icsa_in,file=csa_in,status="old",err=100)
6080 read(icsa_in,*) nconf
6081 read(icsa_in,*) jstart,jend
6082 read(icsa_in,*) nstmax
6083 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6084 read(icsa_in,*) nran0,nran1,irr
6085 read(icsa_in,*) nseed
6086 read(icsa_in,*) ntotal,cut1,cut2
6087 read(icsa_in,*) estop
6088 read(icsa_in,*) icmax,irestart
6089 read(icsa_in,*) ntbankm,dele,difcut
6090 read(icsa_in,*) iref,rmscut,pnccut
6091 read(icsa_in,*) ndiff
6098 end subroutine csa_read
6099 !-----------------------------------------------------------------------------
6100 subroutine initial_write
6103 ! implicit real*8 (a-h,o-z)
6104 ! include 'DIMENSIONS'
6105 ! include 'COMMON.CSA'
6106 ! include 'COMMON.BANK'
6107 ! include 'COMMON.IOUNITS'
6109 ! integer :: ntf,ik,iw_pdb
6110 ! real(kind=8) :: var,ene
6112 open(icsa_seed,file=csa_seed,status="unknown")
6113 write(icsa_seed,*) "seed"
6115 #if defined(AIX) || defined(PGI)
6116 open(icsa_history,file=csa_history,status="unknown",&
6119 open(icsa_history,file=csa_history,status="unknown",&
6122 write(icsa_history,*) nconf
6123 write(icsa_history,*) jstart,jend
6124 write(icsa_history,*) nstmax
6125 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6126 write(icsa_history,*) nran0,nran1,irr
6127 write(icsa_history,*) nseed
6128 write(icsa_history,*) ntotal,cut1,cut2
6129 write(icsa_history,*) estop
6130 write(icsa_history,*) icmax,irestart
6131 write(icsa_history,*) ntbankm,dele,difcut
6132 write(icsa_history,*) iref,rmscut,pnccut
6133 write(icsa_history,*) ndiff
6135 write(icsa_history,*)
6138 open(icsa_bank1,file=csa_bank1,status="unknown")
6139 write(icsa_bank1,*) 0
6143 end subroutine initial_write
6144 !-----------------------------------------------------------------------------
6145 subroutine restart_write
6148 ! implicit real*8 (a-h,o-z)
6149 ! include 'DIMENSIONS'
6150 ! include 'COMMON.IOUNITS'
6151 ! include 'COMMON.CSA'
6152 ! include 'COMMON.BANK'
6154 ! integer :: ntf,ik,iw_pdb
6155 ! real(kind=8) :: var,ene
6157 #if defined(AIX) || defined(PGI)
6158 open(icsa_history,file=csa_history,position="append")
6160 open(icsa_history,file=csa_history,access="append")
6162 write(icsa_history,*)
6163 write(icsa_history,*) "This is restart"
6164 write(icsa_history,*)
6165 write(icsa_history,*) nconf
6166 write(icsa_history,*) jstart,jend
6167 write(icsa_history,*) nstmax
6168 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6169 write(icsa_history,*) nran0,nran1,irr
6170 write(icsa_history,*) nseed
6171 write(icsa_history,*) ntotal,cut1,cut2
6172 write(icsa_history,*) estop
6173 write(icsa_history,*) icmax,irestart
6174 write(icsa_history,*) ntbankm,dele,difcut
6175 write(icsa_history,*) iref,rmscut,pnccut
6176 write(icsa_history,*) ndiff
6177 write(icsa_history,*)
6178 write(icsa_history,*) "irestart is: ", irestart
6180 write(icsa_history,*)
6184 end subroutine restart_write
6185 !-----------------------------------------------------------------------------
6187 !-----------------------------------------------------------------------------
6188 subroutine write_pdb(npdb,titelloc,ee)
6190 ! implicit real*8 (a-h,o-z)
6191 ! include 'DIMENSIONS'
6192 ! include 'COMMON.IOUNITS'
6193 character(len=50) :: titelloc1
6194 character*(*) :: titelloc
6195 character(len=3) :: zahl
6196 character(len=5) :: liczba5
6198 integer :: npdb !,ilen
6202 ! real(kind=8) :: var,ene
6206 if (npdb.lt.1000) then
6207 call numstr(npdb,zahl)
6208 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
6210 if (npdb.lt.10000) then
6211 write(liczba5,'(i1,i4)') 0,npdb
6213 write(liczba5,'(i5)') npdb
6215 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
6217 call pdbout(ee,titelloc1,ipdb)
6220 end subroutine write_pdb
6221 !-----------------------------------------------------------------------------
6223 !-----------------------------------------------------------------------------
6224 subroutine write_thread_summary
6225 ! Thread the sequence through a database of known structures
6226 use control_data, only: refstr
6228 use energy_data, only: n_ene_comp
6230 ! implicit real*8 (a-h,o-z)
6231 ! include 'DIMENSIONS'
6233 use MPI_data !include 'COMMON.INFO'
6236 ! include 'COMMON.CONTROL'
6237 ! include 'COMMON.CHAIN'
6238 ! include 'COMMON.DBASE'
6239 ! include 'COMMON.INTERACT'
6240 ! include 'COMMON.VAR'
6241 ! include 'COMMON.THREAD'
6242 ! include 'COMMON.FFIELD'
6243 ! include 'COMMON.SBRIDGE'
6244 ! include 'COMMON.HEADER'
6245 ! include 'COMMON.NAMES'
6246 ! include 'COMMON.IOUNITS'
6247 ! include 'COMMON.TIME1'
6249 integer,dimension(maxthread) :: ip
6250 real(kind=8),dimension(0:n_ene) :: energia
6252 integer :: i,j,ii,jj,ipj,ik,kk,ist
6253 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
6255 write (iout,'(30x,a/)') &
6256 ' *********** Summary threading statistics ************'
6257 write (iout,'(a)') 'Initial energies:'
6258 write (iout,'(a4,2x,a12,14a14,3a8)') &
6259 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
6260 'RMSnat','NatCONT','NNCONT','RMS'
6261 ! Energy sort patterns
6266 enet=ener(n_ene-1,ip(i))
6269 if (ener(n_ene-1,ip(j)).lt.enet) then
6271 enet=ener(n_ene-1,ip(j))
6283 ist=nres_base(2,ii)+ipatt(2,i)
6285 energia(i)=ener0(kk,i)
6287 etot=ener0(n_ene_comp+1,i)
6288 rmsnat=ener0(n_ene_comp+2,i)
6289 rms=ener0(n_ene_comp+3,i)
6290 frac=ener0(n_ene_comp+4,i)
6291 frac_nn=ener0(n_ene_comp+5,i)
6294 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6295 i,str_nam(ii),ist+1,&
6296 (energia(print_order(kk)),kk=1,nprint_ene),&
6297 etot,rmsnat,frac,frac_nn,rms
6299 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
6300 i,str_nam(ii),ist+1,&
6301 (energia(print_order(kk)),kk=1,nprint_ene),etot
6304 write (iout,'(//a)') 'Final energies:'
6305 write (iout,'(a4,2x,a12,17a14,3a8)') &
6306 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
6307 'RMSnat','NatCONT','NNCONT','RMS'
6311 ist=nres_base(2,ii)+ipatt(2,i)
6313 energia(kk)=ener(kk,ik)
6315 etot=ener(n_ene_comp+1,i)
6316 rmsnat=ener(n_ene_comp+2,i)
6317 rms=ener(n_ene_comp+3,i)
6318 frac=ener(n_ene_comp+4,i)
6319 frac_nn=ener(n_ene_comp+5,i)
6320 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6321 i,str_nam(ii),ist+1,&
6322 (energia(print_order(kk)),kk=1,nprint_ene),&
6323 etot,rmsnat,frac,frac_nn,rms
6325 write (iout,'(/a/)') 'IEXAM array:'
6326 write (iout,'(i5)') nexcl
6328 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6330 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6331 'Max. time for threading step ',max_time_for_thread,&
6332 'Average time for threading step: ',ave_time_for_thread
6334 end subroutine write_thread_summary
6335 !-----------------------------------------------------------------------------
6336 subroutine write_stat_thread(ithread,ipattern,ist)
6338 use energy_data, only: n_ene_comp
6340 ! implicit real*8 (a-h,o-z)
6341 ! include "DIMENSIONS"
6342 ! include "COMMON.CONTROL"
6343 ! include "COMMON.IOUNITS"
6344 ! include "COMMON.THREAD"
6345 ! include "COMMON.FFIELD"
6346 ! include "COMMON.DBASE"
6347 ! include "COMMON.NAMES"
6348 real(kind=8),dimension(0:n_ene) :: energia
6350 integer :: ithread,ipattern,ist,i
6351 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6353 #if defined(AIX) || defined(PGI)
6354 open(istat,file=statname,position='append')
6356 open(istat,file=statname,access='append')
6359 energia(i)=ener(i,ithread)
6361 etot=ener(n_ene_comp+1,ithread)
6362 rmsnat=ener(n_ene_comp+2,ithread)
6363 rms=ener(n_ene_comp+3,ithread)
6364 frac=ener(n_ene_comp+4,ithread)
6365 frac_nn=ener(n_ene_comp+5,ithread)
6366 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6367 ithread,str_nam(ipattern),ist+1,&
6368 (energia(print_order(i)),i=1,nprint_ene),&
6369 etot,rmsnat,frac,frac_nn,rms
6372 end subroutine write_stat_thread
6373 !-----------------------------------------------------------------------------
6375 subroutine readpdb_template(k)
6376 ! Read the PDB file for read_constr_homology with read2sigma
6377 ! and convert the peptide geometry into virtual-chain geometry.
6379 ! include 'DIMENSIONS'
6380 ! include 'COMMON.LOCAL'
6381 ! include 'COMMON.VAR'
6382 ! include 'COMMON.CHAIN'
6383 ! include 'COMMON.INTERACT'
6384 ! include 'COMMON.IOUNITS'
6385 ! include 'COMMON.GEO'
6386 ! include 'COMMON.NAMES'
6387 ! include 'COMMON.CONTROL'
6388 ! include 'COMMON.FRAG'
6389 ! include 'COMMON.SETUP'
6390 use compare_data, only:nhfrag,nbfrag
6391 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6393 logical lprn /.false./,fail
6394 real(kind=8), dimension (3):: e1,e2,e3
6395 real(kind=8) :: dcj,efree_temp
6399 real(kind=8), dimension (3,40) :: sccor
6401 integer, dimension (:), allocatable :: iterter
6402 if(.not.allocated(iterter))allocate(iterter(nres))
6409 write (2,*) "UNRES_PDB",unres_pdb
6417 read (ipdbin,'(a80)',end=10) card
6418 if (card(:3).eq.'END') then
6420 else if (card(:3).eq.'TER') then
6423 itype(ires_old-1,1)=ntyp1
6424 iterter(ires_old-1)=1
6425 #if defined(WHAM_RUN) || defined(CLUSTER)
6426 if (ires_old.lt.nres) then
6428 itype(ires_old,1)=ntyp1
6430 #if defined(WHAM_RUN) || defined(CLUSTER)
6434 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6437 dc(j,ires)=sccor(j,iii)
6440 call sccenter(ires,iii,sccor)
6443 ! Fish out the ATOM cards.
6444 if (index(card(1:4),'ATOM').gt.0) then
6445 read (card(12:16),*) atom
6446 ! write (iout,*) "! ",atom," !",ires
6447 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6448 read (card(23:26),*) ires
6449 read (card(18:20),'(a3)') res
6450 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6451 ! & " ires_old",ires_old
6452 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6453 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6454 if (ires-ishift+ishift1.ne.ires_old) then
6455 ! Calculate the CM of the preceding residue.
6459 dc(j,ires_old)=sccor(j,iii)
6462 call sccenter(ires_old,iii,sccor)
6466 ! Start new residue.
6467 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6470 else if (ibeg.eq.1) then
6471 ! write (iout,*) "BEG ires",ires
6473 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6477 ires=ires-ishift+ishift1
6479 ! write (iout,*) "ishift",ishift," ires",ires,
6480 ! & " ires_old",ires_old
6481 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6483 else if (ibeg.eq.2) then
6485 ishift=-ires_old+ires-1
6487 ! write (iout,*) "New chain started",ires,ishift
6490 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6491 ires=ires-ishift+ishift1
6494 if (res.eq.'ACE' .or. res.eq.'NHE') then
6497 itype(ires,1)=rescode(ires,res,0,1)
6500 ires=ires-ishift+ishift1
6502 ! write (iout,*) "ires_old",ires_old," ires",ires
6503 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6506 ! write (2,*) "ires",ires," res ",res," ity",ity
6507 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6508 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6509 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6510 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6512 write (iout,'(2i3,2x,a,3f8.3)') &
6513 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6517 sccor(j,iii)=c(j,ires)
6519 if (ishift.ne.0) then
6520 ires_ca=ires+ishift-ishift1
6524 ! write (*,*) card(23:27),ires,itype(ires)
6525 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6526 atom.ne.'N' .and. atom.ne.'C' .and.&
6527 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6528 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6529 ! write (iout,*) "sidechain ",atom
6531 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6535 10 if(me.eq.king.or..not.out1file) &
6536 write (iout,'(a,i5)') ' Nres: ',ires
6537 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6540 write(2,*) "tutaj",ires,nres
6542 ! write (iout,*) i,itype(i),itype(i+1)
6543 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6544 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6545 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6546 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6547 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6549 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6550 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6557 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6561 dcj=(c(j,i-2)-c(j,i-3))/2.0
6562 if (dcj.eq.0) dcj=1.23591524223
6567 else !itype(i+1).eq.ntyp1
6569 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6570 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6577 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6581 dcj=(c(j,i+3)-c(j,i+2))/2.0
6582 if (dcj.eq.0) dcj=1.23591524223
6587 endif !itype(i+1).eq.ntyp1
6588 endif !itype.eq.ntyp1
6590 ! Calculate the CM of the last side chain.
6593 dc(j,ires)=sccor(j,iii)
6596 call sccenter(ires,iii,sccor)
6600 if (itype(nres,1).ne.10) then
6604 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6605 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6612 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6616 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6617 if (dcj.eq.0) dcj=1.23591524223
6618 c(j,nres)=c(j,nres-1)+dcj
6619 c(j,2*nres)=c(j,nres)
6630 c(j,2*nres)=c(j,nres)
6632 if (itype(1,1).eq.ntyp1) then
6636 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6637 call refsys(2,3,4,e1,e2,e3,fail)
6644 c(j,1)=c(j,2)-1.9d0*e2(j)
6648 dcj=(c(j,4)-c(j,3))/2.0
6654 ! Copy the coordinates to reference coordinates
6660 ! Calculate internal coordinates.
6661 if (out_template_coord) then
6662 write (iout,'(/a)') &
6663 "Cartesian coordinates of the reference structure"
6664 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6665 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6667 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6668 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6669 (c(j,ires+nres),j=1,3)
6672 ! Calculate internal coordinates.
6673 call int_from_cart(.true.,out_template_coord)
6674 call sc_loc_geom(.false.)
6676 thetaref(i)=theta(i)
6681 dc(j,i)=c(j,i+1)-c(j,i)
6682 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6687 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6688 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6690 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6691 ! & vbld_inv(i+nres)
6696 cref(j,i+nres,1)=c(j,i+nres)
6706 end subroutine readpdb_template
6707 !-----------------------------------------------------------------------------
6709 !-----------------------------------------------------------------------------
6710 end module io_config