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)
834 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
837 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
838 ! dsc(i) = vbldsc0_nucl(1,i)
842 ! dsc_inv(i)=1.0D0/dsc(i)
845 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
846 ! do i=1,ntyp_molec(2)
847 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
848 ! aksc_nucl(j,i),abond0_nucl(j,i),&
849 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
850 ! dsc(i) = vbldsc0(1,i)
854 ! dsc_inv(i)=1.0D0/dsc(i)
859 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
860 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
862 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
864 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
865 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
867 write (iout,'(13x,3f10.5)') &
868 vbldsc0(j,i),aksc(j,i),abond0(j,i)
875 if (.not.allocated(ichargecat)) &
876 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
878 if (oldion.eq.1) then
880 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
881 print *,msc(i,5),restok(i,5)
885 read (iion,*) ncatprotparm
886 allocate(catprm(ncatprotparm,4))
888 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
892 allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
896 read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
899 write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
900 "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
904 write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
905 restyp(i,5),"-",restyp(j,2)
908 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
909 !----------------------------------------------------
910 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
911 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
912 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
913 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
914 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
915 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
925 allocate(liptranene(ntyp))
926 !C reading lipid parameters
927 write (iout,*) "iliptranpar",iliptranpar
929 read(iliptranpar,*) pepliptran
932 read(iliptranpar,*) liptranene(i)
933 print *,liptranene(i)
936 ! water parmaters entalphy
937 allocate(awaterenta(0:400))
938 allocate(bwaterenta(0:400))
939 allocate(cwaterenta(0:400))
940 allocate(dwaterenta(0:400))
941 allocate(awaterentro(0:400))
942 allocate(bwaterentro(0:400))
943 allocate(cwaterentro(0:400))
944 allocate(dwaterentro(0:400))
946 read(iwaterwater,*) mychar
947 read(iwaterwater,*) ract,awaterenta(0),bwaterenta(0),&
948 cwaterenta(0),dwaterenta(0)
950 read(iwaterwater,*) ract,awaterenta(i),bwaterenta(i),&
951 cwaterenta(i),dwaterenta(i)
952 irdiff=int((ract-2.06d0)*50.0d0)+1
953 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
955 ! water parmaters entrophy
956 read(iwaterwater,*) mychar
957 read(iwaterwater,*) ract,awaterentro(0),bwaterentro(0),&
958 cwaterentro(0),dwaterentro(0)
960 read(iwaterwater,*) ract,awaterentro(i),bwaterentro(i),&
961 cwaterentro(i),dwaterentro(i)
962 irdiff=int((ract-2.06d0)*50.0d0)+1
963 if (i.ne.irdiff) print *,"WARTINING",i,ract, irdiff
969 ! Read the parameters of the probability distribution/energy expression
970 ! of the virtual-bond valence angles theta
973 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
974 (bthet(j,i,1,1),j=1,2)
975 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
976 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
977 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
981 athet(1,i,1,-1)=athet(1,i,1,1)
982 athet(2,i,1,-1)=athet(2,i,1,1)
983 bthet(1,i,1,-1)=-bthet(1,i,1,1)
984 bthet(2,i,1,-1)=-bthet(2,i,1,1)
985 athet(1,i,-1,1)=-athet(1,i,1,1)
986 athet(2,i,-1,1)=-athet(2,i,1,1)
987 bthet(1,i,-1,1)=bthet(1,i,1,1)
988 bthet(2,i,-1,1)=bthet(2,i,1,1)
992 athet(1,i,-1,-1)=athet(1,-i,1,1)
993 athet(2,i,-1,-1)=-athet(2,-i,1,1)
994 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
995 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
996 athet(1,i,-1,1)=athet(1,-i,1,1)
997 athet(2,i,-1,1)=-athet(2,-i,1,1)
998 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
999 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1000 athet(1,i,1,-1)=-athet(1,-i,1,1)
1001 athet(2,i,1,-1)=athet(2,-i,1,1)
1002 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1003 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1004 theta0(i)=theta0(-i)
1008 polthet(j,i)=polthet(j,-i)
1011 gthet(j,i)=gthet(j,-i)
1017 if (.not.LaTeX) then
1018 write (iout,'(a)') &
1019 'Parameters of the virtual-bond valence angles:'
1020 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
1021 ' ATHETA0 ',' A1 ',' A2 ',&
1024 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1025 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
1027 write (iout,'(/a/9x,5a/79(1h-))') &
1028 'Parameters of the expression for sigma(theta_c):',&
1029 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
1030 ' ALPH3 ',' SIGMA0C '
1032 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
1033 (polthet(j,i),j=0,3),sigc0(i)
1035 write (iout,'(/a/9x,5a/79(1h-))') &
1036 'Parameters of the second gaussian:',&
1037 ' THETA0 ',' SIGMA0 ',' G1 ',&
1040 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
1041 sig0(i),(gthet(j,i),j=1,3)
1044 write (iout,'(a)') &
1045 'Parameters of the virtual-bond valence angles:'
1046 write (iout,'(/a/9x,5a/79(1h-))') &
1047 'Coefficients of expansion',&
1048 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
1049 ' b1*10^1 ',' b2*10^1 '
1051 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
1052 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
1053 (10*bthet(j,i,1,1),j=1,2)
1055 write (iout,'(/a/9x,5a/79(1h-))') &
1056 'Parameters of the expression for sigma(theta_c):',&
1057 ' alpha0 ',' alph1 ',' alph2 ',&
1058 ' alhp3 ',' sigma0c '
1060 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1061 (polthet(j,i),j=0,3),sigc0(i)
1063 write (iout,'(/a/9x,5a/79(1h-))') &
1064 'Parameters of the second gaussian:',&
1065 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1068 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1069 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1075 ! Read the parameters of Utheta determined from ab initio surfaces
1076 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1078 IF (tor_mode.eq.0) THEN
1079 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1080 ntheterm3,nsingle,ndouble
1081 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1083 !----------------------------------------------------
1084 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1085 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1086 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1087 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1088 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1089 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1090 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1091 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1092 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1093 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1094 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1095 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1096 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1097 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1098 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1099 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1100 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1101 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1102 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1103 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1104 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1105 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1106 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1107 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1109 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1111 ithetyp(i)=-ithetyp(-i)
1114 aa0thet(:,:,:,:)=0.0d0
1115 aathet(:,:,:,:,:)=0.0d0
1116 bbthet(:,:,:,:,:,:)=0.0d0
1117 ccthet(:,:,:,:,:,:)=0.0d0
1118 ddthet(:,:,:,:,:,:)=0.0d0
1119 eethet(:,:,:,:,:,:)=0.0d0
1120 ffthet(:,:,:,:,:,:,:)=0.0d0
1121 ggthet(:,:,:,:,:,:,:)=0.0d0
1123 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1125 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1126 ! VAR:1=non-glicyne non-proline 2=proline
1127 ! VAR:negative values for D-aminoacid
1129 do j=-nthetyp,nthetyp
1130 do k=-nthetyp,nthetyp
1131 read (ithep,'(6a)',end=111,err=111) res1
1132 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1133 ! VAR: aa0thet is variable describing the average value of Foureir
1134 ! VAR: expansion series
1135 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1136 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1137 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1138 read (ithep,*,end=111,err=111) &
1139 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1140 read (ithep,*,end=111,err=111) &
1141 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1142 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1143 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1144 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1146 read (ithep,*,end=111,err=111) &
1147 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1148 ffthet(lll,llll,ll,i,j,k,iblock),&
1149 ggthet(llll,lll,ll,i,j,k,iblock),&
1150 ggthet(lll,llll,ll,i,j,k,iblock),&
1151 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1156 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1157 ! coefficients of theta-and-gamma-dependent terms are zero.
1158 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1159 ! RECOMENTDED AFTER VERSION 3.3)
1163 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1164 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1166 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1167 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1170 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1172 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1175 ! AND COMMENT THE LOOPS BELOW
1179 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1180 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1182 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1183 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1186 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1188 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1193 ! Substitution for D aminoacids from symmetry.
1196 do j=-nthetyp,nthetyp
1197 do k=-nthetyp,nthetyp
1198 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1200 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1204 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1205 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1206 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1207 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1213 ffthet(llll,lll,ll,i,j,k,iblock)= &
1214 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1215 ffthet(lll,llll,ll,i,j,k,iblock)= &
1216 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1217 ggthet(llll,lll,ll,i,j,k,iblock)= &
1218 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1219 ggthet(lll,llll,ll,i,j,k,iblock)= &
1220 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1229 ! Control printout of the coefficients of virtual-bond-angle potentials
1232 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1237 write (iout,'(//4a)') &
1238 'Type ',onelett(i),onelett(j),onelett(k)
1239 write (iout,'(//a,10x,a)') " l","a[l]"
1240 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1241 write (iout,'(i2,1pe15.5)') &
1242 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1244 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1245 "b",l,"c",l,"d",l,"e",l
1247 write (iout,'(i2,4(1pe15.5))') m,&
1248 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1249 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1253 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1254 "f+",l,"f-",l,"g+",l,"g-",l
1257 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1258 ffthet(n,m,l,i,j,k,iblock),&
1259 ffthet(m,n,l,i,j,k,iblock),&
1260 ggthet(n,m,l,i,j,k,iblock),&
1261 ggthet(m,n,l,i,j,k,iblock)
1272 !C here will be the apropriate recalibrating for D-aminoacid
1273 read (ithep,*,end=121,err=121) nthetyp
1274 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1275 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1276 do i=-nthetyp+1,nthetyp-1
1277 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1278 do j=0,nbend_kcc_Tb(i)
1279 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1283 write (iout,'(a)') &
1284 "Parameters of the valence-only potentials"
1285 do i=-nthetyp+1,nthetyp-1
1286 write (iout,'(2a)') "Type ",toronelet(i)
1287 do k=0,nbend_kcc_Tb(i)
1288 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1294 write (2,*) "Start reading THETA_PDB",ithep_pdb
1296 ! write (2,*) 'i=',i
1297 read (ithep_pdb,*,err=111,end=111) &
1298 a0thet(i),(athet(j,i,1,1),j=1,2),&
1299 (bthet(j,i,1,1),j=1,2)
1300 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1301 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1302 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1303 sigc0(i)=sigc0(i)**2
1306 athet(1,i,1,-1)=athet(1,i,1,1)
1307 athet(2,i,1,-1)=athet(2,i,1,1)
1308 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1309 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1310 athet(1,i,-1,1)=-athet(1,i,1,1)
1311 athet(2,i,-1,1)=-athet(2,i,1,1)
1312 bthet(1,i,-1,1)=bthet(1,i,1,1)
1313 bthet(2,i,-1,1)=bthet(2,i,1,1)
1316 a0thet(i)=a0thet(-i)
1317 athet(1,i,-1,-1)=athet(1,-i,1,1)
1318 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1319 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1320 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1321 athet(1,i,-1,1)=athet(1,-i,1,1)
1322 athet(2,i,-1,1)=-athet(2,-i,1,1)
1323 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1324 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1325 athet(1,i,1,-1)=-athet(1,-i,1,1)
1326 athet(2,i,1,-1)=athet(2,-i,1,1)
1327 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1328 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1329 theta0(i)=theta0(-i)
1333 polthet(j,i)=polthet(j,-i)
1336 gthet(j,i)=gthet(j,-i)
1339 write (2,*) "End reading THETA_PDB"
1343 !--------------- Reading theta parameters for nucleic acid-------
1344 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1345 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1346 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1347 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1348 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1349 nthetyp_nucl+1,nthetyp_nucl+1))
1350 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1351 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1352 nthetyp_nucl+1,nthetyp_nucl+1))
1353 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1354 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1355 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1356 nthetyp_nucl+1,nthetyp_nucl+1))
1357 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1358 nthetyp_nucl+1,nthetyp_nucl+1))
1359 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1360 nthetyp_nucl+1,nthetyp_nucl+1))
1361 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1362 nthetyp_nucl+1,nthetyp_nucl+1))
1363 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1364 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1365 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1366 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1367 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1368 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1370 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1371 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1373 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1375 aa0thet_nucl(:,:,:)=0.0d0
1376 aathet_nucl(:,:,:,:)=0.0d0
1377 bbthet_nucl(:,:,:,:,:)=0.0d0
1378 ccthet_nucl(:,:,:,:,:)=0.0d0
1379 ddthet_nucl(:,:,:,:,:)=0.0d0
1380 eethet_nucl(:,:,:,:,:)=0.0d0
1381 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1382 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1387 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1388 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1389 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1390 read (ithep_nucl,*,end=111,err=111) &
1391 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1392 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1393 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1394 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1395 read (ithep_nucl,*,end=111,err=111) &
1396 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1397 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1398 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1403 !-------------------------------------------
1404 allocate(nlob(ntyp1)) !(ntyp1)
1405 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1406 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1407 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1414 gaussc(:,:,:,:)=0.0D0
1418 ! Read the parameters of the probability distribution/energy expression
1419 ! of the side chains.
1422 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1426 dsc_inv(i)=1.0D0/dsc(i)
1437 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1438 ((blower(k,l,1),l=1,k),k=1,3)
1439 censc(1,1,-i)=censc(1,1,i)
1440 censc(2,1,-i)=censc(2,1,i)
1441 censc(3,1,-i)=-censc(3,1,i)
1443 read (irotam,*,end=112,err=112) bsc(j,i)
1444 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1445 ((blower(k,l,j),l=1,k),k=1,3)
1446 censc(1,j,-i)=censc(1,j,i)
1447 censc(2,j,-i)=censc(2,j,i)
1448 censc(3,j,-i)=-censc(3,j,i)
1449 ! BSC is amplitude of Gaussian
1456 akl=akl+blower(k,m,j)*blower(l,m,j)
1460 if (((k.eq.3).and.(l.ne.3)) &
1461 .or.((l.eq.3).and.(k.ne.3))) then
1462 gaussc(k,l,j,-i)=-akl
1463 gaussc(l,k,j,-i)=-akl
1465 gaussc(k,l,j,-i)=akl
1466 gaussc(l,k,j,-i)=akl
1475 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1478 if (nlobi.gt.0) then
1480 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1481 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1482 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1483 'log h',(bsc(j,i),j=1,nlobi)
1484 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1485 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1487 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1488 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1491 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1492 write (iout,'(a,f10.4,4(16x,f10.4))') &
1493 'Center ',(bsc(j,i),j=1,nlobi)
1494 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1503 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1504 ! added by Urszula Kozlowska 07/11/2007
1506 !el Maximum number of SC local term fitting function coefficiants
1507 !el integer,parameter :: maxsccoef=65
1509 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1512 read (irotam,*,end=112,err=112)
1514 read (irotam,*,end=112,err=112)
1517 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1522 allocate(nterm_scend(2,ntyp))
1523 allocate(arotam_end(0:6,2,ntyp))
1526 read (irotam_end,*) ijunk
1527 !c write (iout,*) "ijunk",ijunk
1531 read (irotam_end,'(a)')
1532 read (irotam_end,*) nterm_scend(j,i)
1533 !c write (iout,*) "i",i," j",j," nterm",nterm_scend(j,i)
1534 do k=0,nterm_scend(j,i)
1535 read (irotam_end,*) ijunk,arotam_end(k,j,i)
1536 !c write (iout,*) "k",k," arotam",arotam_end(k,j,i)
1542 write (iout,'(a)') &
1543 "Parameters of the local potentials of sidechain ends"
1545 write (iout,'(5x,9x,2hp-,a3,6x,9x,a3,2h-p)')&
1546 restyp(i,1),restyp(i,1)
1547 do j=0,max0(nterm_scend(1,i),nterm_scend(2,i))
1548 write (iout,'(i5,2f20.10)') &
1549 j,arotam_end(j,1,i),arotam_end(j,2,i)
1556 !---------reading nucleic acid parameters for rotamers-------------------
1557 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1558 do i=1,ntyp_molec(2)
1559 read (irotam_nucl,*,end=112,err=112)
1561 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1567 write (iout,*) "Base rotamer parameters"
1568 do i=1,ntyp_molec(2)
1569 write (iout,'(a)') restyp(i,2)
1570 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1575 ! Read the parameters of the probability distribution/energy expression
1576 ! of the side chains.
1578 write (2,*) "Start reading ROTAM_PDB"
1580 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1584 dsc_inv(i)=1.0D0/dsc(i)
1595 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1596 ((blower(k,l,1),l=1,k),k=1,3)
1598 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1599 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1600 ((blower(k,l,j),l=1,k),k=1,3)
1607 akl=akl+blower(k,m,j)*blower(l,m,j)
1617 write (2,*) "End reading ROTAM_PDB"
1623 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1624 !C interaction energy of the Gly, Ala, and Pro prototypes.
1626 read (ifourier,*) nloctyp
1627 SPLIT_FOURIERTOR = nloctyp.lt.0
1628 nloctyp = iabs(nloctyp)
1629 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1630 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1631 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1632 !C allocate(ctilde(2,2,nres))
1633 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1634 !C allocate(gtb1(2,nres))
1635 !C allocate(gtb2(2,nres))
1636 !C allocate(cc(2,2,nres))
1637 !C allocate(dd(2,2,nres))
1638 !C allocate(ee(2,2,nres))
1639 !C allocate(gtcc(2,2,nres))
1640 !C allocate(gtdd(2,2,nres))
1641 !C allocate(gtee(2,2,nres))
1644 allocate(itype2loc(-ntyp1:ntyp1))
1645 allocate(iloctyp(-nloctyp:nloctyp))
1646 allocate(bnew1(3,2,-nloctyp:nloctyp))
1647 allocate(bnew2(3,2,-nloctyp:nloctyp))
1648 allocate(ccnew(3,2,-nloctyp:nloctyp))
1649 allocate(ddnew(3,2,-nloctyp:nloctyp))
1650 allocate(e0new(3,-nloctyp:nloctyp))
1651 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1652 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1653 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1654 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1655 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1656 allocate(e0newtor(3,-nloctyp:nloctyp))
1657 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1670 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1671 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1672 itype2loc(ntyp1)=nloctyp
1673 iloctyp(nloctyp)=ntyp1
1675 itype2loc(-i)=-itype2loc(i)
1678 allocate(iloctyp(-nloctyp:nloctyp))
1679 allocate(itype2loc(-ntyp1:ntyp1))
1686 iloctyp(-i)=-iloctyp(i)
1688 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1689 !c write (iout,*) "nloctyp",nloctyp,
1690 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1691 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1692 !c write (iout,*) "nloctyp",nloctyp,
1693 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1696 !c write (iout,*) "NEWCORR",i
1697 read (ifourier,*,end=115,err=115)
1700 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1703 !c write (iout,*) "NEWCORR BNEW1"
1704 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1707 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1710 !c write (iout,*) "NEWCORR BNEW2"
1711 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1713 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1714 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1716 !c write (iout,*) "NEWCORR CCNEW"
1717 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1719 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1720 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1722 !c write (iout,*) "NEWCORR DDNEW"
1723 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1727 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1731 !c write (iout,*) "NEWCORR EENEW1"
1732 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1734 read (ifourier,*,end=115,err=115) e0new(ii,i)
1736 !c write (iout,*) (e0new(ii,i),ii=1,3)
1738 !c write (iout,*) "NEWCORR EENEW"
1741 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1742 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1743 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1744 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1749 bnew1(ii,1,-i)= bnew1(ii,1,i)
1750 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1751 bnew2(ii,1,-i)= bnew2(ii,1,i)
1752 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1755 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1756 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1757 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1758 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1759 ccnew(ii,1,-i)=ccnew(ii,1,i)
1760 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1761 ddnew(ii,1,-i)=ddnew(ii,1,i)
1762 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1764 e0new(1,-i)= e0new(1,i)
1765 e0new(2,-i)=-e0new(2,i)
1766 e0new(3,-i)=-e0new(3,i)
1768 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1769 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1770 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1771 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1775 write (iout,'(a)') "Coefficients of the multibody terms"
1776 do i=-nloctyp+1,nloctyp-1
1777 write (iout,*) "Type: ",onelet(iloctyp(i))
1778 write (iout,*) "Coefficients of the expansion of B1"
1780 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1782 write (iout,*) "Coefficients of the expansion of B2"
1784 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1786 write (iout,*) "Coefficients of the expansion of C"
1787 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1788 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1789 write (iout,*) "Coefficients of the expansion of D"
1790 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1791 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1792 write (iout,*) "Coefficients of the expansion of E"
1793 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1796 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1801 IF (SPLIT_FOURIERTOR) THEN
1803 !c write (iout,*) "NEWCORR TOR",i
1804 read (ifourier,*,end=115,err=115)
1807 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1810 !c write (iout,*) "NEWCORR BNEW1 TOR"
1811 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1814 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1817 !c write (iout,*) "NEWCORR BNEW2 TOR"
1818 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1820 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1821 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1823 !c write (iout,*) "NEWCORR CCNEW TOR"
1824 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1826 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1827 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1829 !c write (iout,*) "NEWCORR DDNEW TOR"
1830 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1834 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1838 !c write (iout,*) "NEWCORR EENEW1 TOR"
1839 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1841 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1843 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1845 !c write (iout,*) "NEWCORR EENEW TOR"
1848 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1849 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1850 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1851 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1856 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1857 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1858 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1859 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1862 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1863 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1864 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1865 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1867 e0newtor(1,-i)= e0newtor(1,i)
1868 e0newtor(2,-i)=-e0newtor(2,i)
1869 e0newtor(3,-i)=-e0newtor(3,i)
1871 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1872 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1873 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1874 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1878 write (iout,'(a)') &
1879 "Single-body coefficients of the torsional potentials"
1880 do i=-nloctyp+1,nloctyp-1
1881 write (iout,*) "Type: ",onelet(iloctyp(i))
1882 write (iout,*) "Coefficients of the expansion of B1tor"
1884 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1885 j,(bnew1tor(k,j,i),k=1,3)
1887 write (iout,*) "Coefficients of the expansion of B2tor"
1889 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1890 j,(bnew2tor(k,j,i),k=1,3)
1892 write (iout,*) "Coefficients of the expansion of Ctor"
1893 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1894 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1895 write (iout,*) "Coefficients of the expansion of Dtor"
1896 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1897 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1898 write (iout,*) "Coefficients of the expansion of Etor"
1899 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1902 write (iout,'(1hE,2i1,2f10.5)') &
1903 j,k,(eenewtor(l,j,k,i),l=1,2)
1909 do i=-nloctyp+1,nloctyp-1
1912 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1917 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1921 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1922 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1923 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1924 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1927 ENDIF !SPLIT_FOURIER_TOR
1929 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1930 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1931 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1932 allocate(b(13,-nloctyp-1:nloctyp+1))
1934 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1936 read (ifourier,*,end=115,err=115)
1937 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1939 write (iout,*) 'Type ',onelet(iloctyp(i))
1940 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1954 !c B1tilde(1,i) = b(3)
1955 !c! B1tilde(2,i) =-b(5)
1956 !c! B1tilde(1,-i) =-b(3)
1957 !c! B1tilde(2,-i) =b(5)
1958 !c! b1tilde(1,i)=0.0d0
1959 !c b1tilde(2,i)=0.0d0
1964 !cc B1tilde(1,i) = b(3,i)
1965 !cc B1tilde(2,i) =-b(5,i)
1966 !c B1tilde(1,-i) =-b(3,i)
1967 !c B1tilde(2,-i) =b(5,i)
1968 !cc b1tilde(1,i)=0.0d0
1969 !cc b1tilde(2,i)=0.0d0
1970 !cc B2(1,i) = b(2,i)
1971 !cc B2(2,i) = b(4,i)
1973 !c B2(2,-i) =-b(4,i)
1977 CCold(1,1,i)= b(7,i)
1978 CCold(2,2,i)=-b(7,i)
1979 CCold(2,1,i)= b(9,i)
1980 CCold(1,2,i)= b(9,i)
1981 CCold(1,1,-i)= b(7,i)
1982 CCold(2,2,-i)=-b(7,i)
1983 CCold(2,1,-i)=-b(9,i)
1984 CCold(1,2,-i)=-b(9,i)
1989 !c Ctilde(1,1,i)= CCold(1,1,i)
1990 !c Ctilde(1,2,i)= CCold(1,2,i)
1991 !c Ctilde(2,1,i)=-CCold(2,1,i)
1992 !c Ctilde(2,2,i)=-CCold(2,2,i)
1997 !c Ctilde(1,1,i)= CCold(1,1,i)
1998 !c Ctilde(1,2,i)= CCold(1,2,i)
1999 !c Ctilde(2,1,i)=-CCold(2,1,i)
2000 !c Ctilde(2,2,i)=-CCold(2,2,i)
2002 !c Ctilde(1,1,i)=0.0d0
2003 !c Ctilde(1,2,i)=0.0d0
2004 !c Ctilde(2,1,i)=0.0d0
2005 !c Ctilde(2,2,i)=0.0d0
2006 DDold(1,1,i)= b(6,i)
2007 DDold(2,2,i)=-b(6,i)
2008 DDold(2,1,i)= b(8,i)
2009 DDold(1,2,i)= b(8,i)
2010 DDold(1,1,-i)= b(6,i)
2011 DDold(2,2,-i)=-b(6,i)
2012 DDold(2,1,-i)=-b(8,i)
2013 DDold(1,2,-i)=-b(8,i)
2018 !c Dtilde(1,1,i)= DD(1,1,i)
2019 !c Dtilde(1,2,i)= DD(1,2,i)
2020 !c Dtilde(2,1,i)=-DD(2,1,i)
2021 !c Dtilde(2,2,i)=-DD(2,2,i)
2023 !c Dtilde(1,1,i)=0.0d0
2024 !c Dtilde(1,2,i)=0.0d0
2025 !c Dtilde(2,1,i)=0.0d0
2026 !c Dtilde(2,2,i)=0.0d0
2027 EEold(1,1,i)= b(10,i)+b(11,i)
2028 EEold(2,2,i)=-b(10,i)+b(11,i)
2029 EEold(2,1,i)= b(12,i)-b(13,i)
2030 EEold(1,2,i)= b(12,i)+b(13,i)
2031 EEold(1,1,-i)= b(10,i)+b(11,i)
2032 EEold(2,2,-i)=-b(10,i)+b(11,i)
2033 EEold(2,1,-i)=-b(12,i)+b(13,i)
2034 EEold(1,2,-i)=-b(12,i)-b(13,i)
2035 write(iout,*) "TU DOCHODZE"
2041 !c ee(2,1,i)=ee(1,2,i)
2046 "Coefficients of the cumulants (independent of valence angles)"
2047 do i=-nloctyp+1,nloctyp-1
2048 write (iout,*) 'Type ',onelet(iloctyp(i))
2050 write(iout,'(2f10.5)') B(3,i),B(5,i)
2052 write(iout,'(2f10.5)') B(2,i),B(4,i)
2055 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
2059 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
2063 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
2072 ! Read torsional parameters in old format
2074 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
2076 read (itorp,*,end=113,err=113) ntortyp,nterm_old
2077 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
2078 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2080 !el from energy module--------
2081 allocate(v1(nterm_old,ntortyp,ntortyp))
2082 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
2083 !el---------------------------
2088 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
2094 write (iout,'(/a/)') 'Torsional constants:'
2097 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
2098 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2104 ! Read torsional parameters
2106 IF (TOR_MODE.eq.0) THEN
2107 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2108 read (itorp,*,end=113,err=113) ntortyp
2109 !el from energy module---------
2110 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2111 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2113 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2114 allocate(vlor2(maxlor,ntortyp,ntortyp))
2115 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2116 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2118 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2119 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2120 !el---------------------------
2123 !el---------------------------
2125 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2127 itortyp(i)=-itortyp(-i)
2129 itortyp(ntyp1)=ntortyp
2130 itortyp(-ntyp1)=-ntortyp
2132 write (iout,*) 'ntortyp',ntortyp
2134 do j=-ntortyp+1,ntortyp-1
2135 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2137 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2138 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2141 do k=1,nterm(i,j,iblock)
2142 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2144 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2145 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2146 v0ij=v0ij+si*v1(k,i,j,iblock)
2148 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2149 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2150 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2152 do k=1,nlor(i,j,iblock)
2153 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2154 vlor2(k,i,j),vlor3(k,i,j)
2155 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2158 v0(-i,-j,iblock)=v0ij
2164 write (iout,'(/a/)') 'Torsional constants:'
2166 do i=-ntortyp,ntortyp
2167 do j=-ntortyp,ntortyp
2168 write (iout,*) 'ityp',i,' jtyp',j
2169 write (iout,*) 'Fourier constants'
2170 do k=1,nterm(i,j,iblock)
2171 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2174 write (iout,*) 'Lorenz constants'
2175 do k=1,nlor(i,j,iblock)
2176 write (iout,'(3(1pe15.5))') &
2177 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2183 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2185 ! 6/23/01 Read parameters for double torsionals
2187 !el from energy module------------
2188 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2189 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2190 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2191 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2192 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2193 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2194 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2195 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2196 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2197 !---------------------------------
2201 do j=-ntortyp+1,ntortyp-1
2202 do k=-ntortyp+1,ntortyp-1
2203 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2204 ! write (iout,*) "OK onelett",
2207 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2208 .or. t3.ne.toronelet(k)) then
2209 write (iout,*) "Error in double torsional parameter file",&
2212 call MPI_Finalize(Ierror)
2214 stop "Error in double torsional parameter file"
2216 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2217 ntermd_2(i,j,k,iblock)
2218 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2219 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2220 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2221 ntermd_1(i,j,k,iblock))
2222 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2223 ntermd_1(i,j,k,iblock))
2224 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2225 ntermd_1(i,j,k,iblock))
2226 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2227 ntermd_1(i,j,k,iblock))
2228 ! Martix of D parameters for one dimesional foureir series
2229 do l=1,ntermd_1(i,j,k,iblock)
2230 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2231 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2232 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2233 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2234 ! write(iout,*) "whcodze" ,
2235 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2237 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2238 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2239 v2s(m,l,i,j,k,iblock),&
2240 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2241 ! Martix of D parameters for two dimesional fourier series
2242 do l=1,ntermd_2(i,j,k,iblock)
2244 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2245 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2246 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2247 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2256 write (iout,*) 'Constants for double torsionals'
2259 do j=-ntortyp+1,ntortyp-1
2260 do k=-ntortyp+1,ntortyp-1
2261 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2262 ' nsingle',ntermd_1(i,j,k,iblock),&
2263 ' ndouble',ntermd_2(i,j,k,iblock)
2265 write (iout,*) 'Single angles:'
2266 do l=1,ntermd_1(i,j,k,iblock)
2267 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2268 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2269 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2270 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2273 write (iout,*) 'Pairs of angles:'
2274 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2275 do l=1,ntermd_2(i,j,k,iblock)
2276 write (iout,'(i5,20f10.5)') &
2277 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2280 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2281 do l=1,ntermd_2(i,j,k,iblock)
2282 write (iout,'(i5,20f10.5)') &
2283 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2284 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2294 itype2loc(i)=itortyp(i)
2298 ELSE IF (TOR_MODE.eq.1) THEN
2300 !C read valence-torsional parameters
2301 read (itorp,*,end=121,err=121) ntortyp
2303 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2305 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2306 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2309 itype2loc(i)=itortyp(i)
2313 itortyp(i)=-itortyp(-i)
2315 do i=-ntortyp+1,ntortyp-1
2316 do j=-ntortyp+1,ntortyp-1
2317 !C first we read the cos and sin gamma parameters
2318 read (itorp,'(13x,a)',end=121,err=121) string
2319 write (iout,*) i,j,string
2320 read (itorp,*,end=121,err=121) &
2321 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2322 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2323 do k=1,nterm_kcc(j,i)
2324 do l=1,nterm_kcc_Tb(j,i)
2325 do ll=1,nterm_kcc_Tb(j,i)
2326 read (itorp,*,end=121,err=121) ii,jj,kk, &
2327 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2335 !c AL 4/8/16: Calculate coefficients from one-body parameters
2337 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2338 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2339 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2340 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2341 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2344 print *,i,itortyp(i)
2345 itortyp(i)=itype2loc(i)
2348 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2350 do i=-ntortyp+1,ntortyp-1
2351 do j=-ntortyp+1,ntortyp-1
2354 do k=1,nterm_kcc_Tb(j,i)
2355 do l=1,nterm_kcc_Tb(j,i)
2356 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2357 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2358 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2359 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2362 do k=1,nterm_kcc_Tb(j,i)
2363 do l=1,nterm_kcc_Tb(j,i)
2365 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2366 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2367 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2368 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2370 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2371 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2372 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2373 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2377 !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)
2381 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2386 if (tor_mode.gt.0 .and. lprint) then
2387 !c Print valence-torsional parameters
2388 write (iout,'(a)') &
2389 "Parameters of the valence-torsional potentials"
2390 do i=-ntortyp+1,ntortyp-1
2391 do j=-ntortyp+1,ntortyp-1
2392 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2393 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2394 do k=1,nterm_kcc(j,i)
2395 do l=1,nterm_kcc_Tb(j,i)
2396 do ll=1,nterm_kcc_Tb(j,i)
2397 write (iout,'(3i5,2f15.4)')&
2398 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2406 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2407 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2408 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2409 !el from energy module---------
2410 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2411 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2413 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2414 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2415 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2416 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2418 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2419 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2420 !el---------------------------
2423 !el--------------------
2424 read (itorp_nucl,*,end=113,err=113) &
2425 (itortyp_nucl(i),i=1,ntyp_molec(2))
2426 ! print *,itortyp_nucl(:)
2427 !c write (iout,*) 'ntortyp',ntortyp
2430 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2431 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2434 do k=1,nterm_nucl(i,j)
2435 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2436 v0ij=v0ij+si*v1_nucl(k,i,j)
2439 do k=1,nlor_nucl(i,j)
2440 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2441 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2442 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2448 ! Read of Side-chain backbone correlation parameters
2449 ! Modified 11 May 2012 by Adasko
2452 read (isccor,*,end=119,err=119) nsccortyp
2454 !el from module energy-------------
2455 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2456 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2457 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2458 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2459 !-----------------------------------
2461 !el from module energy-------------
2462 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2464 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2466 isccortyp(i)=-isccortyp(-i)
2468 iscprol=isccortyp(20)
2469 ! write (iout,*) 'ntortyp',ntortyp
2471 !c maxinter is maximum interaction sites
2472 !el from module energy---------
2473 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2474 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2475 -nsccortyp:nsccortyp))
2476 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2477 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2478 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2479 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2480 !-----------------------------------
2484 read (isccor,*,end=119,err=119) &
2485 nterm_sccor(i,j),nlor_sccor(i,j)
2491 nterm_sccor(-i,j)=nterm_sccor(i,j)
2492 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2493 nterm_sccor(i,-j)=nterm_sccor(i,j)
2494 do k=1,nterm_sccor(i,j)
2495 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2497 if (j.eq.iscprol) then
2498 if (i.eq.isccortyp(10)) then
2499 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2500 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2502 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2503 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2504 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2505 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2506 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2507 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2508 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2509 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2512 if (i.eq.isccortyp(10)) then
2513 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2514 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2516 if (j.eq.isccortyp(10)) then
2517 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2518 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2520 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2521 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2522 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2523 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2524 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2525 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2529 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2530 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2531 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2532 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2535 do k=1,nlor_sccor(i,j)
2536 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2537 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2538 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2539 (1+vlor3sccor(k,i,j)**2)
2541 v0sccor(l,i,j)=v0ijsccor
2542 v0sccor(l,-i,j)=v0ijsccor1
2543 v0sccor(l,i,-j)=v0ijsccor2
2544 v0sccor(l,-i,-j)=v0ijsccor3
2550 !el from module energy-------------
2551 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2553 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2554 ! write (iout,*) 'ntortyp',ntortyp
2556 !c maxinter is maximum interaction sites
2557 !el from module energy---------
2558 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2559 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2560 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2561 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2562 !-----------------------------------
2566 read (isccor,*,end=119,err=119) &
2567 nterm_sccor(i,j),nlor_sccor(i,j)
2571 do k=1,nterm_sccor(i,j)
2572 print *,"test",i,j,k,l
2573 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2575 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2578 do k=1,nlor_sccor(i,j)
2579 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2580 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2581 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2582 (1+vlor3sccor(k,i,j)**2)
2584 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2593 write (iout,'(/a/)') 'Torsional constants:'
2596 write (iout,*) 'ityp',i,' jtyp',j
2597 write (iout,*) 'Fourier constants'
2598 do k=1,nterm_sccor(i,j)
2599 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2601 write (iout,*) 'Lorenz constants'
2602 do k=1,nlor_sccor(i,j)
2603 write (iout,'(3(1pe15.5))') &
2604 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2611 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2612 ! interaction energy of the Gly, Ala, and Pro prototypes.
2615 ! Read electrostatic-interaction parameters
2620 write (iout,'(/a)') 'Electrostatic interaction constants:'
2621 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2622 'IT','JT','APP','BPP','AEL6','AEL3'
2624 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2625 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2626 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2627 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2632 app (i,j)=epp(i,j)*rri*rri
2633 bpp (i,j)=-2.0D0*epp(i,j)*rri
2634 ael6(i,j)=elpp6(i,j)*4.2D0**6
2635 ael3(i,j)=elpp3(i,j)*4.2D0**3
2637 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2643 ! Read side-chain interaction parameters.
2645 !el from module energy - COMMON.INTERACT-------
2646 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2647 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2648 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2650 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2651 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2652 allocate(epslip(ntyp,ntyp))
2660 !--------------------------------
2662 read (isidep,*,end=117,err=117) ipot,expon
2663 if (ipot.lt.1 .or. ipot.gt.5) then
2664 write (iout,'(2a)') 'Error while reading SC interaction',&
2665 'potential file - unknown potential type.'
2667 call MPI_Finalize(Ierror)
2673 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2674 ', exponents are ',expon,2*expon
2675 ! goto (10,20,30,30,40) ipot
2677 !----------------------- LJ potential ---------------------------------
2679 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2680 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2681 (sigma0(i),i=1,ntyp)
2683 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2684 write (iout,'(a/)') 'The epsilon array:'
2685 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2686 write (iout,'(/a)') 'One-body parameters:'
2687 write (iout,'(a,4x,a)') 'residue','sigma'
2688 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2691 !----------------------- LJK potential --------------------------------
2693 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2694 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2695 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2697 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2698 write (iout,'(a/)') 'The epsilon array:'
2699 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2700 write (iout,'(/a)') 'One-body parameters:'
2701 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2702 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2706 !---------------------- GB or BP potential -----------------------------
2709 ! print *,"I AM in SCELE",scelemode
2710 if (scelemode.eq.0) then
2712 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2714 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2715 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2716 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2717 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2719 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2722 ! For the GB potential convert sigma'**2 into chi'
2725 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2729 write (iout,'(/a/)') 'Parameters of the BP potential:'
2730 write (iout,'(a/)') 'The epsilon array:'
2731 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2732 write (iout,'(/a)') 'One-body parameters:'
2733 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2735 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2736 chip(i),alp(i),i=1,ntyp)
2739 ! print *,ntyp,"NTYP"
2740 allocate(icharge(ntyp1))
2741 ! print *,ntyp,icharge(i)
2743 read (isidep,*) (icharge(i),i=1,ntyp)
2744 print *,ntyp,icharge(i)
2745 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2746 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2747 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2748 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2749 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2750 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2751 allocate(epsintab(ntyp,ntyp))
2752 allocate(dtail(2,ntyp,ntyp))
2753 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2754 allocate(wqdip(2,ntyp,ntyp))
2755 allocate(wstate(4,ntyp,ntyp))
2756 allocate(dhead(2,2,ntyp,ntyp))
2757 allocate(nstate(ntyp,ntyp))
2758 allocate(debaykap(ntyp,ntyp))
2760 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2761 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2765 ! write (*,*) "Im in ALAB", i, " ", j
2767 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2768 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2769 chis(i,j),chis(j,i), & !2 w tej linii
2770 nstate(i,j),(wstate(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
2771 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),& ! 4 w tej linii
2772 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2773 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2774 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), & ! 5 w tej linii
2775 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2776 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2777 IF ((LaTeX).and.(i.gt.24)) then
2778 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2779 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), & !6 w tej linii
2780 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), & !6 w tej linii
2781 chis(i,j),chis(j,i) !2 w tej linii
2783 print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i), wqdip(1,i,j)
2787 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2791 IF ((LaTeX).and.(i.gt.24)) then
2792 write (2,'(2a4,1h&,14(f8.2,1h&),23(f8.2,1h&))') restyp(i,1),restyp(j,1), &
2793 dhead(1,1,i,j),dhead(2,1,i,j),& ! 2 w tej linii
2794 dtail(1,i,j),dtail(2,i,j), & ! 2 w tej lini
2795 epshead(i,j),sig0head(i,j), & ! 2 w tej linii
2796 rborn(i,j),rborn(j,i), & ! 3 w tej linii
2797 alphapol(i,j),alphapol(j,i), & ! 2 w tej linii
2798 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j) !8 w tej linii
2805 sigma(i,j) = sigma(j,i)
2806 sigmap1(i,j)=sigmap1(j,i)
2807 sigmap2(i,j)=sigmap2(j,i)
2808 sigiso1(i,j)=sigiso1(j,i)
2809 sigiso2(i,j)=sigiso2(j,i)
2810 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2811 nstate(i,j) = nstate(j,i)
2812 dtail(1,i,j) = dtail(2,j,i)
2813 dtail(2,i,j) = dtail(1,j,i)
2815 alphasur(k,i,j) = alphasur(k,j,i)
2816 wstate(k,i,j) = wstate(k,j,i)
2817 alphiso(k,i,j) = alphiso(k,j,i)
2820 dhead(2,1,i,j) = dhead(1,1,j,i)
2821 dhead(2,2,i,j) = dhead(1,2,j,i)
2822 dhead(1,1,i,j) = dhead(2,1,j,i)
2823 dhead(1,2,i,j) = dhead(2,2,j,i)
2825 epshead(i,j) = epshead(j,i)
2826 sig0head(i,j) = sig0head(j,i)
2829 wqdip(k,i,j) = wqdip(k,j,i)
2832 wquad(i,j) = wquad(j,i)
2833 epsintab(i,j) = epsintab(j,i)
2834 debaykap(i,j)=debaykap(j,i)
2835 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2842 !--------------------- GBV potential -----------------------------------
2844 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2845 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2846 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2847 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2849 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2850 write (iout,'(a/)') 'The epsilon array:'
2851 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2852 write (iout,'(/a)') 'One-body parameters:'
2853 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2854 's||/s_|_^2',' chip ',' alph '
2855 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2856 sigii(i),chip(i),alp(i),i=1,ntyp)
2859 write(iout,*)"Wrong ipot"
2865 !-----------------------------------------------------------------------
2866 ! Calculate the "working" parameters of SC interactions.
2868 !el from module energy - COMMON.INTERACT-------
2869 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2870 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2871 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2872 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2873 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2874 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2876 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2882 if (scelemode.eq.0) then
2891 sc_aa_tube_par(:)=0.0d0
2892 sc_bb_tube_par(:)=0.0d0
2894 !--------------------------------
2899 epslip(i,j)=epslip(j,i)
2902 if (scelemode.eq.0) then
2905 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2906 sigma(j,i)=sigma(i,j)
2907 rs0(i,j)=dwa16*sigma(i,j)
2912 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2913 'Working parameters of the SC interactions:',&
2914 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2919 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2921 ! print *,"SIGMA ZLA?",sigma(i,j)
2929 sigeps=dsign(1.0D0,epsij)
2931 aa_aq(i,j)=epsij*rrij*rrij
2932 ! print *,"ADASKO",epsij,rrij,expon
2933 bb_aq(i,j)=-sigeps*epsij*rrij
2934 aa_aq(j,i)=aa_aq(i,j)
2935 bb_aq(j,i)=bb_aq(i,j)
2936 epsijlip=epslip(i,j)
2937 sigeps=dsign(1.0D0,epsijlip)
2938 epsijlip=dabs(epsijlip)
2939 aa_lip(i,j)=epsijlip*rrij*rrij
2940 bb_lip(i,j)=-sigeps*epsijlip*rrij
2941 aa_lip(j,i)=aa_lip(i,j)
2942 bb_lip(j,i)=bb_lip(i,j)
2943 !C write(iout,*) aa_lip
2944 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2945 sigt1sq=sigma0(i)**2
2946 sigt2sq=sigma0(j)**2
2949 ratsig1=sigt2sq/sigt1sq
2950 ratsig2=1.0D0/ratsig1
2951 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2952 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2953 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2957 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2958 sigmaii(i,j)=rsum_max
2959 sigmaii(j,i)=rsum_max
2961 ! sigmaii(i,j)=r0(i,j)
2962 ! sigmaii(j,i)=r0(i,j)
2964 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2965 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2966 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2967 augm(i,j)=epsij*r_augm**(2*expon)
2968 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2975 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2976 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2977 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2982 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2983 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2984 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2985 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2986 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2987 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2988 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2989 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2990 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2991 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2992 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2993 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2994 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2995 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2996 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2997 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
3006 read (isidep_nucl,*) ipot_nucl
3007 ! print *,"TU?!",ipot_nucl
3008 if (ipot_nucl.eq.1) then
3009 do i=1,ntyp_molec(2)
3010 do j=i,ntyp_molec(2)
3011 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
3012 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
3016 do i=1,ntyp_molec(2)
3017 do j=i,ntyp_molec(2)
3018 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
3019 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
3020 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
3024 ! rpp(1,1)=2**(1.0/6.0)*5.16158
3025 do i=1,ntyp_molec(2)
3026 do j=i,ntyp_molec(2)
3027 rrij=sigma_nucl(i,j)
3031 epsij=4*eps_nucl(i,j)
3032 sigeps=dsign(1.0D0,epsij)
3034 aa_nucl(i,j)=epsij*rrij*rrij
3035 bb_nucl(i,j)=-sigeps*epsij*rrij
3036 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
3037 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
3038 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
3039 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
3040 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
3041 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
3044 aa_nucl(i,j)=aa_nucl(j,i)
3045 bb_nucl(i,j)=bb_nucl(j,i)
3046 ael3_nucl(i,j)=ael3_nucl(j,i)
3047 ael6_nucl(i,j)=ael6_nucl(j,i)
3048 ael63_nucl(i,j)=ael63_nucl(j,i)
3049 ael32_nucl(i,j)=ael32_nucl(j,i)
3050 elpp3_nucl(i,j)=elpp3_nucl(j,i)
3051 elpp6_nucl(i,j)=elpp6_nucl(j,i)
3052 elpp63_nucl(i,j)=elpp63_nucl(j,i)
3053 elpp32_nucl(i,j)=elpp32_nucl(j,i)
3054 eps_nucl(i,j)=eps_nucl(j,i)
3055 sigma_nucl(i,j)=sigma_nucl(j,i)
3056 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
3060 write(iout,*) "tube param"
3061 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
3062 ccavtubpep,dcavtubpep,tubetranenepep
3063 sigmapeptube=sigmapeptube**6
3064 sigeps=dsign(1.0D0,epspeptube)
3065 epspeptube=dabs(epspeptube)
3066 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
3067 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
3068 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
3070 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
3071 ccavtub(i),dcavtub(i),tubetranene(i)
3072 sigmasctube=sigmasctube**6
3073 sigeps=dsign(1.0D0,epssctube)
3074 epssctube=dabs(epssctube)
3075 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
3076 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
3077 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
3079 !-----------------READING SC BASE POTENTIALS-----------------------------
3080 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
3081 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
3082 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
3083 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
3084 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
3085 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
3086 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
3087 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
3088 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
3089 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
3090 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
3091 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
3092 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
3093 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
3094 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
3095 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
3097 write (iout,*) "ESCBASEPARM"
3098 do i=1,ntyp_molec(1)
3099 do j=1,ntyp_molec(2) ! without U then we will take T for U
3100 ! write (*,*) "Im in ", i, " ", j
3101 read(isidep_scbase,*) &
3102 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3103 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3104 ! write(*,*) "eps",eps_scbase(i,j)
3105 read(isidep_scbase,*) &
3106 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3107 chis_scbase(i,j,1),chis_scbase(i,j,2)
3108 read(isidep_scbase,*) &
3109 dhead_scbasei(i,j), &
3110 dhead_scbasej(i,j), &
3111 rborn_scbasei(i,j),rborn_scbasej(i,j)
3112 read(isidep_scbase,*) &
3113 (wdipdip_scbase(k,i,j),k=1,3), &
3114 (wqdip_scbase(k,i,j),k=1,2)
3115 read(isidep_scbase,*) &
3116 alphapol_scbase(i,j), &
3117 epsintab_scbase(i,j)
3118 if (chi_scbase(i,j,2).gt.0.9) chi_scbase(i,j,2)=0.9
3119 if (chi_scbase(i,j,1).gt.0.9) chi_scbase(i,j,1)=0.9
3120 if (chipp_scbase(i,j,2).gt.0.9) chipp_scbase(i,j,2)=0.9
3121 if (chipp_scbase(i,j,1).gt.0.9) chipp_scbase(i,j,1)=0.9
3122 if (chi_scbase(i,j,2).lt.-0.9) chi_scbase(i,j,2)=-0.9
3123 if (chi_scbase(i,j,1).lt.-0.9) chi_scbase(i,j,1)=-0.9
3124 if (chipp_scbase(i,j,2).lt.-0.9) chipp_scbase(i,j,2)=-0.9
3125 if (chipp_scbase(i,j,1).lt.-0.9) chipp_scbase(i,j,1)=-0.9
3127 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3128 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3129 write(*,*) "eps",eps_scbase(i,j)
3131 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3132 chis_scbase(i,j,1),chis_scbase(i,j,2)
3134 dhead_scbasei(i,j), &
3135 dhead_scbasej(i,j), &
3136 rborn_scbasei(i,j),rborn_scbasej(i,j)
3138 (wdipdip_scbase(k,i,j),k=1,3), &
3139 (wqdip_scbase(k,i,j),k=1,2)
3141 alphapol_scbase(i,j), &
3142 epsintab_scbase(i,j)
3147 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
3148 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
3149 write(*,*) "eps",eps_scbase(i,j)
3151 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
3152 chis_scbase(i,j,1),chis_scbase(i,j,2)
3154 dhead_scbasei(i,j), &
3155 dhead_scbasej(i,j), &
3156 rborn_scbasei(i,j),rborn_scbasej(i,j)
3158 (wdipdip_scbase(k,i,j),k=1,3), &
3159 (wqdip_scbase(k,i,j),k=1,2)
3161 alphapol_scbase(i,j), &
3162 epsintab_scbase(i,j)
3165 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
3166 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
3168 do i=1,ntyp_molec(1)
3169 do j=1,ntyp_molec(2)-1
3170 epsij=eps_scbase(i,j)
3171 rrij=sigma_scbase(i,j)
3176 sigeps=dsign(1.0D0,epsij)
3178 aa_scbase(i,j)=epsij*rrij*rrij
3179 bb_scbase(i,j)=-sigeps*epsij*rrij
3182 !-----------------READING PEP BASE POTENTIALS-------------------
3183 allocate(eps_pepbase(ntyp_molec(2)))
3184 allocate(sigma_pepbase(ntyp_molec(2)))
3185 allocate(chi_pepbase(ntyp_molec(2),2))
3186 allocate(chipp_pepbase(ntyp_molec(2),2))
3187 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3188 allocate(sigmap1_pepbase(ntyp_molec(2)))
3189 allocate(sigmap2_pepbase(ntyp_molec(2)))
3190 allocate(chis_pepbase(ntyp_molec(2),2))
3191 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3193 write (iout,*) "EPEPBASEPARM"
3195 do j=1,ntyp_molec(2) ! without U then we will take T for U
3196 write (*,*) "Im in ", i, " ", j
3197 read(isidep_pepbase,*) &
3198 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3199 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3200 if (chi_pepbase(j,2).gt.0.9) chi_pepbase(j,2)=0.9
3201 if (chi_pepbase(j,1).gt.0.9) chi_pepbase(j,1)=0.9
3202 if (chipp_pepbase(j,2).gt.0.9) chipp_pepbase(j,2)=0.9
3203 if (chipp_pepbase(j,1).gt.0.9) chipp_pepbase(j,1)=0.9
3204 if (chi_pepbase(j,2).lt.-0.9) chi_pepbase(j,2)=-0.9
3205 if (chi_pepbase(j,1).lt.-0.9) chi_pepbase(j,1)=-0.9
3206 if (chipp_pepbase(j,2).lt.-0.9) chipp_pepbase(j,2)=-0.9
3207 if (chipp_pepbase(j,1).lt.-0.9) chipp_pepbase(j,1)=-0.9
3209 write(*,*) "eps",eps_pepbase(j)
3210 read(isidep_pepbase,*) &
3211 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3212 chis_pepbase(j,1),chis_pepbase(j,2)
3213 read(isidep_pepbase,*) &
3214 (wdipdip_pepbase(k,j),k=1,3)
3215 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3216 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3218 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3219 chis_pepbase(j,1),chis_pepbase(j,2)
3221 (wdipdip_pepbase(k,j),k=1,3)
3225 write(iout,*) eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3226 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3228 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3229 chis_pepbase(j,1),chis_pepbase(j,2)
3231 (wdipdip_pepbase(k,j),k=1,3)
3233 allocate(aa_pepbase(ntyp_molec(2)))
3234 allocate(bb_pepbase(ntyp_molec(2)))
3236 do j=1,ntyp_molec(2)-1
3237 epsij=eps_pepbase(j)
3238 rrij=sigma_pepbase(j)
3243 sigeps=dsign(1.0D0,epsij)
3245 aa_pepbase(j)=epsij*rrij*rrij
3246 bb_pepbase(j)=-sigeps*epsij*rrij
3248 !--------------READING SC PHOSPHATE-------------------------------------
3249 allocate(eps_scpho(ntyp_molec(1)))
3250 allocate(sigma_scpho(ntyp_molec(1)))
3251 allocate(chi_scpho(ntyp_molec(1),2))
3252 allocate(chipp_scpho(ntyp_molec(1),2))
3253 allocate(alphasur_scpho(4,ntyp_molec(1)))
3254 allocate(sigmap1_scpho(ntyp_molec(1)))
3255 allocate(sigmap2_scpho(ntyp_molec(1)))
3256 allocate(chis_scpho(ntyp_molec(1),2))
3257 allocate(wqq_scpho(ntyp_molec(1)))
3258 allocate(wqdip_scpho(2,ntyp_molec(1)))
3259 allocate(alphapol_scpho(ntyp_molec(1)))
3260 allocate(epsintab_scpho(ntyp_molec(1)))
3261 allocate(dhead_scphoi(ntyp_molec(1)))
3262 allocate(rborn_scphoi(ntyp_molec(1)))
3263 allocate(rborn_scphoj(ntyp_molec(1)))
3264 allocate(alphi_scpho(ntyp_molec(1)))
3268 do j=1,ntyp_molec(1) ! without U then we will take T for U
3269 write (*,*) "Im in scpho ", i, " ", j
3270 read(isidep_scpho,*) &
3271 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3272 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3273 write(*,*) "eps",eps_scpho(j)
3274 read(isidep_scpho,*) &
3275 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3276 chis_scpho(j,1),chis_scpho(j,2)
3277 read(isidep_scpho,*) &
3278 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3279 read(isidep_scpho,*) &
3280 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3282 if (chi_scpho(j,2).gt.0.9) chi_scpho(j,2)=0.9
3283 if (chi_scpho(j,1).gt.0.9) chi_scpho(j,1)=0.9
3284 if (chipp_scpho(j,2).gt.0.9) chipp_scpho(j,2)=0.9
3285 if (chipp_scpho(j,1).gt.0.9) chipp_scpho(j,1)=0.9
3286 if (chi_scpho(j,2).lt.-0.9) chi_scpho(j,2)=-0.9
3287 if (chi_scpho(j,1).lt.-0.9) chi_scpho(j,1)=-0.9
3288 if (chipp_scpho(j,2).lt.-0.9) chipp_scpho(j,2)=-0.9
3289 if (chipp_scpho(j,1).lt.-0.9) chipp_scpho(j,1)=-0.9
3293 allocate(aa_scpho(ntyp_molec(1)))
3294 allocate(bb_scpho(ntyp_molec(1)))
3296 do j=1,ntyp_molec(1)
3303 sigeps=dsign(1.0D0,epsij)
3305 aa_scpho(j)=epsij*rrij*rrij
3306 bb_scpho(j)=-sigeps*epsij*rrij
3310 read(isidep_peppho,*) &
3311 eps_peppho,sigma_peppho
3312 read(isidep_peppho,*) &
3313 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3314 read(isidep_peppho,*) &
3315 (wqdip_peppho(k),k=1,2)
3323 sigeps=dsign(1.0D0,epsij)
3325 aa_peppho=epsij*rrij*rrij
3326 bb_peppho=-sigeps*epsij*rrij
3329 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3334 ! Define the SC-p interaction constants (hard-coded; old style)
3337 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3339 ! aad(i,1)=0.3D0*4.0D0**12
3340 ! Following line for constants currently implemented
3341 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3342 aad(i,1)=1.5D0*4.0D0**12
3343 ! aad(i,1)=0.17D0*5.6D0**12
3345 ! "Soft" SC-p repulsion
3347 ! Following line for constants currently implemented
3348 ! aad(i,1)=0.3D0*4.0D0**6
3349 ! "Hard" SC-p repulsion
3350 bad(i,1)=3.0D0*4.0D0**6
3351 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3360 ! 8/9/01 Read the SC-p interaction constants from file
3363 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3366 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3367 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3368 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3369 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3373 write (iout,*) "Parameters of SC-p interactions:"
3375 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3376 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3381 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3383 do i=1,ntyp_molec(2)
3384 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3386 do i=1,ntyp_molec(2)
3387 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3388 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3390 r0pp=1.12246204830937298142*5.16158
3396 ! Define the constants of the disulfide bridge
3400 ! Old arbitrary potential - commented out.
3405 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3406 ! energy surface of diethyl disulfide.
3407 ! A. Liwo and U. Kozlowska, 11/24/03
3423 allocate(ichargelipid(ntyp_molec(4)))
3424 allocate(lip_angle_force(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3425 allocate(lip_angle_angle(ntyp_molec(4),ntyp_molec(4),ntyp_molec(4)))
3426 allocate(lip_bond(ntyp_molec(4),ntyp_molec(4)))
3427 allocate(lip_eps(ntyp_molec(4),ntyp_molec(4)))
3428 allocate(lip_sig(ntyp_molec(4),ntyp_molec(4)))
3429 kjtokcal=0.2390057361
3431 !HERE THE MASS of MARTINI
3432 write(*,*) "before MARTINI PARAM"
3433 do i=1,ntyp_molec(4)
3439 msc(ntyp_molec(4)+1,4)=0.1d0
3440 !relative dielectric constant = 15 for implicit screening
3441 k_coulomb_lip=332.0d0/15.0d0
3442 !kbond = 1250 kJ/(mol*nm*2)
3443 kbondlip=1250.0d0*kjtokcal/100.0d0
3445 lip_angle_force=0.0d0
3446 if (DRY_MARTINI.gt.0) then
3447 lip_angle_force(3,12,12)=35.0*kjtokcal!*krad2
3448 lip_angle_force(3,12,18)=35.0*kjtokcal!*krad2
3449 lip_angle_force(3,18,16)=35.0*kjtokcal!*krad2
3450 lip_angle_force(12,18,16)=35.0*kjtokcal!*krad2
3451 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3452 lip_angle_force(16,18,18)=35.0*kjtokcal!*krad2
3453 lip_angle_force(12,18,18)=35.0*kjtokcal!*krad2
3454 lip_angle_force(18,18,18)=35.0*kjtokcal!*krad2
3456 lip_angle_force(3,12,12)=25.0*kjtokcal!*krad2
3457 lip_angle_force(3,12,18)=25.0*kjtokcal!*krad2
3458 lip_angle_force(3,18,16)=25.0*kjtokcal!*krad2
3459 lip_angle_force(12,18,16)=25.0*kjtokcal!*krad2
3460 lip_angle_force(18,16,18)=45.0*kjtokcal!*krad2
3461 lip_angle_force(16,18,18)=25.0*kjtokcal!*krad2
3462 lip_angle_force(12,18,18)=25.0*kjtokcal!*krad2
3463 lip_angle_force(18,18,18)=25.0*kjtokcal!*krad2
3465 lip_angle_angle=0.0d0
3466 lip_angle_angle(3,12,12)=120.0/krad
3467 lip_angle_angle(3,12,18)=180.0/krad
3468 lip_angle_angle(3,18,16)=180.0/krad
3469 lip_angle_angle(12,18,16)=180.0/krad
3470 lip_angle_angle(18,16,18)=120.0/krad
3471 lip_angle_angle(16,18,18)=180.0/krad
3472 lip_angle_angle(12,18,18)=180.0/krad
3473 lip_angle_angle(18,18,18)=180.0/krad
3474 read(ilipbond,*) temp1
3476 read(ilipbond,*) temp1, lip_bond(i,1), &
3477 lip_bond(i,2),lip_bond(i,3),lip_bond(i,4),lip_bond(i,5), &
3478 lip_bond(i,6),lip_bond(i,7),lip_bond(i,8),lip_bond(i,9), &
3479 lip_bond(i,10),lip_bond(i,11),lip_bond(i,12),lip_bond(i,13), &
3480 lip_bond(i,14),lip_bond(i,15),lip_bond(i,16),lip_bond(i,17), &
3483 lip_bond(i,j)=lip_bond(i,j)*10
3487 read(ilipnonbond,*) (ichargelipid(i),i=1,ntyp_molec(4))
3488 read(ilipnonbond,*) temp1
3490 read(ilipnonbond,*) temp1, lip_eps(i,1), &
3491 lip_eps(i,2),lip_eps(i,3),lip_eps(i,4),lip_eps(i,5), &
3492 lip_eps(i,6),lip_eps(i,7),lip_eps(i,8),lip_eps(i,9), &
3493 lip_eps(i,10),lip_eps(i,11),lip_eps(i,12),lip_eps(i,13), &
3494 lip_eps(i,14),lip_eps(i,15),lip_eps(i,16),lip_eps(i,17), &
3496 ! write(*,*) i, lip_eps(i,18)
3498 lip_eps(i,j)=lip_eps(i,j)*kjtokcal
3501 read(ilipnonbond,*) temp1
3503 read(ilipnonbond,*) temp1,lip_sig(i,1), &
3504 lip_sig(i,2),lip_sig(i,3),lip_sig(i,4),lip_sig(i,5), &
3505 lip_sig(i,6),lip_sig(i,7),lip_sig(i,8),lip_sig(i,9), &
3506 lip_sig(i,10),lip_sig(i,11),lip_sig(i,12),lip_sig(i,13), &
3507 lip_sig(i,14),lip_sig(i,15),lip_sig(i,16),lip_sig(i,17), &
3510 lip_sig(i,j)=lip_sig(i,j)*10.0
3513 write(*,*) "after MARTINI PARAM"
3517 allocate(alphapolcat(ntyp,-1:ntyp_molec(5)),epsheadcat(ntyp,-1:ntyp_molec(5)),sig0headcat(ntyp,-1:ntyp_molec(5)))
3518 allocate(alphapolcat2(ntyp,-1:ntyp_molec(5)))
3519 allocate(sigiso1cat(ntyp,-1:ntyp_molec(5)),rborn1cat(ntyp,-1:ntyp_molec(5)),rborn2cat(ntyp,-1:ntyp_molec(5)),sigmap1cat(ntyp,-1:ntyp_molec(5)))
3520 allocate(sigmap2cat(ntyp,-1:ntyp_molec(5)),sigiso2cat(ntyp,-1:ntyp_molec(5)))
3521 allocate(chis1cat(ntyp,-1:ntyp_molec(5)),chis2cat(ntyp,-1:ntyp_molec(5)),wquadcat(ntyp,-1:ntyp_molec(5)),chipp1cat(ntyp,-1:ntyp_molec(5)),chipp2cat(ntyp,-1:ntyp_molec(5)))
3522 allocate(epsintabcat(ntyp,-1:ntyp_molec(5)))
3523 allocate(dtailcat(2,ntyp,-1:ntyp_molec(5)))
3524 allocate(alphasurcat(4,ntyp,-1:ntyp_molec(5)),alphisocat(4,ntyp,-1:ntyp_molec(5)))
3525 allocate(wqdipcat(2,ntyp,-1:ntyp_molec(5)))
3526 allocate(wstatecat(4,ntyp,-1:ntyp_molec(5)))
3527 allocate(dheadcat(2,2,ntyp,-1:ntyp_molec(5)))
3528 allocate(nstatecat(ntyp,-1:ntyp_molec(5)))
3529 allocate(debaykapcat(ntyp,-1:ntyp_molec(5)))
3531 if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,-1:ntyp1))
3532 if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,-1:ntyp1))
3533 ! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
3534 if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
3535 if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,-1:ntyp1)) !(ntyp,ntyp)
3538 if (.not.allocated(ichargecat))&
3539 allocate (ichargecat(-ntyp_molec(5):ntyp_molec(5)))
3540 write(*,*) "before ions",oldion
3543 ! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
3544 if (oldion.eq.0) then
3545 if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
3546 allocate(icharge(1:ntyp1))
3547 read(iion,*) (icharge(i),i=1,ntyp)
3551 print *,ntyp_molec(5)
3552 do i=-ntyp_molec(5),ntyp_molec(5)
3553 read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
3554 print *,msc(i,5),restok(i,5)
3560 do j=-1,ntyp_molec(5)-1 ! this is without Zn will be modified for ALL tranistion metals
3563 ! do j=1,ntyp_molec(5)
3564 ! write (*,*) "Im in ALAB", i, " ", j
3566 epscat(i,j),sigmacat(i,j), &
3567 ! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
3568 chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), & !6
3570 (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&!12
3571 ! chiscat(i,j),chiscat(j,i), &
3572 chis1cat(i,j),chis2cat(i,j), &
3574 nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !19 !5 w tej lini - 1 integer pierwszy
3575 dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&!23
3576 dtailcat(1,i,j),dtailcat(2,i,j), &
3577 epsheadcat(i,j),sig0headcat(i,j), &!27
3579 ! rborncat(i,j),rborncat(j,i),&
3580 rborn1cat(i,j),rborn2cat(i,j),&
3581 (wqdipcat(k,i,j),k=1,2), &!31
3582 alphapolcat(i,j),alphapolcat2(j,i), &!33
3583 (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
3585 if (chi1cat(i,j).gt.0.9) write (*,*) "WTF ANISO", i,j, chi1cat(i,j)
3586 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
3588 ! write (iout,*) 'i= ', i, ' j= ', j
3589 ! write (iout,*) 'epsi0= ', epscat(1,j)
3590 ! write (iout,*) 'sigma0= ', sigmacat(1,j)
3591 ! write (iout,*) 'chi(i,j)= ', chicat(1,j)
3592 ! write (iout,*) 'chi(j,i)= ', chicat(j,1)
3593 ! write (iout,*) 'chip(i,j)= ', chippcat(1,j)
3594 ! write (iout,*) 'chip(j,i)= ', chippcat(j,1)
3595 ! write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
3596 ! write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
3597 ! write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
3598 ! write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
3599 ! write (iout,*) 'sig1= ', sigmap1cat(1,j)
3600 ! write (iout,*) 'chis(i,j)= ', chiscat(1,j)
3601 ! write (iout,*) 'chis(j,i)= ', chiscat(j,1)
3602 ! write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
3603 ! write (iout,*) 'a1= ', rborncat(j,1)
3604 ! write (iout,*) 'a2= ', rborncat(1,j)
3605 ! write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
3606 ! write (iout,*) 'alphapol1= ', alphapolcat(1,j)
3607 ! write (iout,*) 'w1= ', wqdipcat(1,1,j)
3608 ! write (iout,*) 'w2= ', wqdipcat(2,1,j)
3612 ! If ((i.eq.1).and.(j.eq.27)) then
3613 ! write (iout,*) 'SEP'
3614 ! Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3615 ! Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3620 allocate(aa_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)),&
3621 bb_aq_cat(-ntyp:ntyp,-1:ntyp_molec(5)))
3623 do j=-1,ntyp_molec(5)
3628 sigeps=dsign(1.0D0,epsij)
3630 aa_aq_cat(i,j)=epsij*rrij*rrij
3631 bb_aq_cat(i,j)=-sigeps*epsij*rrij
3636 do j=1,ntyp_molec(5)-1
3638 write (iout,*) 'i= ', i, ' j= ', j
3639 write (iout,*) 'epsi0= ', epscat(i,j)
3640 write (iout,*) 'sigma0= ', sigmacat(i,j)
3641 write (iout,*) 'chi1= ', chi1cat(i,j)
3642 write (iout,*) 'chi1= ', chi2cat(i,j)
3643 write (iout,*) 'chip1= ', chipp1cat(i,j)
3644 write (iout,*) 'chip2= ', chipp2cat(i,j)
3645 write (iout,*) 'alphasur1= ', alphasurcat(1,i,j)
3646 write (iout,*) 'alphasur2= ', alphasurcat(2,i,j)
3647 write (iout,*) 'alphasur3= ', alphasurcat(3,i,j)
3648 write (iout,*) 'alphasur4= ', alphasurcat(4,i,j)
3649 write (iout,*) 'sig1= ', sigmap1cat(i,j)
3650 write (iout,*) 'sig2= ', sigmap2cat(i,j)
3651 write (iout,*) 'chis1= ', chis1cat(i,j)
3652 write (iout,*) 'chis1= ', chis2cat(i,j)
3653 write (iout,*) 'nstatecat(i,j)= ', nstatecat(i,j)
3654 write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,i,j)
3655 write (iout,*) 'dhead= ', dheadcat(1,1,i,j)
3656 write (iout,*) 'dhead2= ', dheadcat(1,2,i,j)
3657 write (iout,*) 'a1= ', rborn1cat(i,j)
3658 write (iout,*) 'a2= ', rborn2cat(i,j)
3659 write (iout,*) 'epsin= ', epsintabcat(i,j), epsintabcat(j,i)
3660 write (iout,*) 'alphapol1= ', alphapolcat(i,j)
3661 write (iout,*) 'alphapol2= ', alphapolcat2(i,j)
3662 write (iout,*) 'w1= ', wqdipcat(1,i,j)
3663 write (iout,*) 'w2= ', wqdipcat(2,i,j)
3664 write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(i,j)
3667 If ((i.eq.1).and.(j.eq.27)) then
3668 write (iout,*) 'SEP'
3669 Write (iout,*) 'w1= ', wqdipcat(1,1,27)
3670 Write (iout,*) 'w2= ', wqdipcat(2,1,27)
3677 ! read number of Zn2+
3678 ! here two denotes the Zn2+ and Cu2+
3679 write(iout,*) "before TRANPARM"
3680 allocate(aomicattr(0:3,2))
3681 allocate(athetacattran(0:6,5,2))
3682 allocate(agamacattran(3,5,2))
3683 allocate(acatshiftdsc(5,2))
3684 allocate(bcatshiftdsc(5,2))
3685 allocate(demorsecat(5,2))
3686 allocate(alphamorsecat(5,2))
3687 allocate(x0catleft(5,2))
3688 allocate(x0catright(5,2))
3689 allocate(x0cattrans(5,2))
3690 allocate(ntrantyp(2))
3691 do i=1,1 ! currently only Zn2+
3693 read(iiontran,*) ntrantyp(i)
3696 !ASP| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3697 !CYS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0left x0right x0transi
3698 !GLU| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -0.5
3699 !HIS| a1 a2 a3 aa0 aa1 aa2 aa3 aa4 aa5 aa6 ad bd De alpha x0 -1 -.5
3700 read(iiontran,*) (aomicattr(j,i),j=0,3)
3702 read (iiontran,*,err=123,end=123) (agamacattran(k,j,i),k=1,3),&
3703 (athetacattran(k,j,i),k=0,6),acatshiftdsc(j,i),bcatshiftdsc(j,i),&
3704 demorsecat(j,i),alphamorsecat(j,i),x0catleft(j,i),x0catright(j,i),&
3709 write (iout,'(/a)') "Disulfide bridge parameters:"
3710 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3711 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3712 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3713 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3717 !------------MARTINI-PROTEIN-parameters-------------------------
3718 allocate(alphapolmart(ntyp,ntyp),epsheadmart(ntyp,ntyp_molec(4)),sig0headmart(ntyp,ntyp_molec(4)))
3719 allocate(alphapolmart2(ntyp,ntyp))
3720 allocate(sigiso1mart(ntyp,ntyp_molec(4)),rborn1mart(ntyp,ntyp_molec(4)),rborn2mart(ntyp,ntyp_molec(4)),sigmap1mart(ntyp,ntyp_molec(4)))
3721 allocate(sigmap2mart(ntyp,ntyp_molec(4)),sigiso2mart(ntyp,ntyp_molec(4)))
3722 allocate(chis1mart(ntyp,ntyp_molec(4)),chis2mart(ntyp,ntyp_molec(4)),wquadmart(ntyp,ntyp_molec(4)),chipp1mart(ntyp,ntyp_molec(4)),chipp2mart(ntyp,ntyp_molec(4)))
3723 allocate(epsintabmart(ntyp,ntyp_molec(4)))
3724 allocate(dtailmart(2,ntyp,ntyp_molec(4)))
3725 allocate(alphasurmart(4,ntyp,ntyp_molec(4)),alphisomart(4,ntyp,ntyp_molec(4)))
3726 allocate(wqdipmart(2,ntyp,ntyp_molec(4)))
3727 allocate(wstatemart(4,ntyp,ntyp_molec(4)))
3728 allocate(dheadmart(2,2,ntyp,ntyp_molec(4)))
3729 allocate(nstatemart(ntyp,ntyp_molec(4)))
3730 allocate(debaykapmart(ntyp,ntyp_molec(4)))
3732 if (.not.allocated(epsmart)) allocate (epsmart(0:ntyp1,ntyp1))
3733 if (.not.allocated(sigmamart)) allocate(sigmamart(0:ntyp1,ntyp1))
3734 ! if (.not.allocated(chimart)) allomarte(chimart(ntyp1,ntyp1)) !(ntyp,ntyp)
3735 if (.not.allocated(chi1mart)) allocate(chi1mart(ntyp1,ntyp1)) !(ntyp,ntyp)
3736 if (.not.allocated(chi2mart)) allocate(chi2mart(ntyp1,ntyp1)) !(ntyp,ntyp)
3739 do i=1,ntyp-3 ! there are phosporylated missing
3740 do j=1,ntyp_molec(4) ! this is without Zn will be modified for ALL tranistion metals
3741 ! do j=1,ntyp_molec(5)
3742 print *,"lipmart",i,j
3743 ! write (*,*) "Im in ALAB", i, " ", j
3745 epsmart(i,j),sigmamart(i,j), &
3746 ! chimart(i,j),chimart(j,i),chippmart(i,j),chippmart(j,i), &
3747 chi1mart(i,j),chi2mart(i,j),chipp1mart(i,j),chipp2mart(i,j), & !6
3749 (alphasurmart(k,i,j),k=1,4),sigmap1mart(i,j),sigmap2mart(i,j),&!12
3750 ! chismart(i,j),chismart(j,i), &
3751 chis1mart(i,j),chis2mart(i,j), &
3753 nstatemart(i,j),(wstatemart(k,i,j),k=1,4), & !19 !5 w tej lini - 1 integer pierwszy
3754 dheadmart(1,1,i,j),dheadmart(1,2,i,j),dheadmart(2,1,i,j),dheadmart(2,2,i,j),&!23
3755 dtailmart(1,i,j),dtailmart(2,i,j), &
3756 epsheadmart(i,j),sig0headmart(i,j), &!27
3758 ! rbornmart(i,j),rbornmart(j,i),&
3759 rborn1mart(i,j),rborn2mart(i,j),&
3760 (wqdipmart(k,i,j),k=1,2), &!31
3761 alphapolmart(i,j),alphapolmart2(j,i), &!33
3762 (alphisomart(k,i,j),k=1,4),sigiso1mart(i,j),sigiso2mart(i,j),epsintabmart(i,j),debaykapmart(i,j)
3765 allocate(aa_aq_mart(-ntyp:ntyp,ntyp_molec(4)),&
3766 bb_aq_mart(-ntyp:ntyp,ntyp_molec(4)))
3767 do i=1,ntyp-3 ! still no phophorylated residues
3768 do j=1,ntyp_molec(4)
3773 sigeps=dsign(1.0D0,epsij)
3775 aa_aq_mart(i,j)=epsij*rrij*rrij
3776 bb_aq_mart(i,j)=-sigeps*epsij*rrij
3790 if (shield_mode.gt.0) then
3791 pi=4.0D0*datan(1.0D0)
3792 !C VSolvSphere the volume of solving sphere
3794 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3795 !C there will be no distinction between proline peptide group and normal peptide
3796 !C group in case of shielding parameters
3797 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3798 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3799 write (iout,*) VSolvSphere,VSolvSphere_div
3800 !C long axis of side chain
3802 long_r_sidechain(i)=vbldsc0(1,i)
3803 ! if (scelemode.eq.0) then
3804 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3805 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3807 ! short_r_sidechain(i)=sigma(i,i)
3809 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3816 111 write (iout,*) "Error reading bending energy parameters."
3818 112 write (iout,*) "Error reading rotamer energy parameters."
3820 113 write (iout,*) "Error reading torsional energy parameters."
3822 114 write (iout,*) "Error reading double torsional energy parameters."
3824 115 write (iout,*) &
3825 "Error reading cumulant (multibody energy) parameters."
3827 116 write (iout,*) "Error reading electrostatic energy parameters."
3829 117 write (iout,*) "Error reading side chain interaction parameters."
3831 118 write (iout,*) "Error reading SCp interaction parameters."
3833 119 write (iout,*) "Error reading SCCOR parameters"
3835 121 write (iout,*) "Error in Czybyshev parameters"
3837 123 write(iout,*) "Error in transition metal parameters"
3840 call MPI_Finalize(Ierror)
3844 end subroutine parmread
3846 !-----------------------------------------------------------------------------
3848 !-----------------------------------------------------------------------------
3849 subroutine printmat(ldim,m,n,iout,key,a)
3852 character(len=3),dimension(n) :: key
3853 real(kind=8),dimension(ldim,n) :: a
3855 integer :: i,j,k,m,iout,nlim
3859 write (iout,1000) (key(k),k=i,nlim)
3861 1000 format (/5x,8(6x,a3))
3862 1020 format (/80(1h-)/)
3864 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3867 1010 format (a3,2x,8(f9.4))
3869 end subroutine printmat
3870 !-----------------------------------------------------------------------------
3872 !-----------------------------------------------------------------------------
3874 ! Read the PDB file and convert the peptide geometry into virtual-chain
3877 use energy_data, only: itype,istype
3881 ! use control, only: rescode,sugarcode
3882 ! implicit real*8 (a-h,o-z)
3883 ! include 'DIMENSIONS'
3884 ! include 'COMMON.LOCAL'
3885 ! include 'COMMON.VAR'
3886 ! include 'COMMON.CHAIN'
3887 ! include 'COMMON.INTERACT'
3888 ! include 'COMMON.IOUNITS'
3889 ! include 'COMMON.GEO'
3890 ! include 'COMMON.NAMES'
3891 ! include 'COMMON.CONTROL'
3892 ! include 'COMMON.DISTFIT'
3893 ! include 'COMMON.SETUP'
3894 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3896 logical :: lprn=.true.,fail
3897 real(kind=8),dimension(3) :: e1,e2,e3
3898 real(kind=8) :: dcj,efree_temp
3899 character(len=3) :: seq,res,res2
3900 character(len=5) :: atom
3901 character(len=80) :: card
3902 real(kind=8),dimension(3,40) :: sccor
3903 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3904 integer :: isugar,molecprev,firstion
3905 character*1 :: sugar
3907 real(kind=8),dimension(3) :: ccc
3909 integer,dimension(2,maxres/3) :: hfrag_alloc
3910 integer,dimension(4,maxres/3) :: bfrag_alloc
3911 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3912 real(kind=8),dimension(:,:), allocatable :: c_temporary
3913 integer,dimension(:,:) , allocatable :: itype_temporary
3914 integer,dimension(:),allocatable :: istype_temp
3921 ! write (2,*) "UNRES_PDB",unres_pdb
3941 !-----------------------------
3942 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3943 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3944 if(.not. allocated(istype)) allocate(istype(maxres))
3946 read (ipdbin,'(a80)',end=10) card
3947 write (iout,'(a)') card
3948 if (card(:5).eq.'HELIX') then
3951 read(card(22:25),*) hfrag(1,nhfrag)
3952 read(card(34:37),*) hfrag(2,nhfrag)
3954 if (card(:5).eq.'SHEET') then
3957 read(card(24:26),*) bfrag(1,nbfrag)
3958 read(card(35:37),*) bfrag(2,nbfrag)
3959 !rc----------------------------------------
3960 !rc to be corrected !!!
3961 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3962 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3963 !rc----------------------------------------
3965 if (card(:3).eq.'END') then
3967 else if (card(:3).eq.'TER') then
3972 itype(ires_old,molecule)=ntyp1_molec(molecule)
3973 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3974 nres_molec(molecule)=nres_molec(molecule)+2
3975 ! if (molecule.eq.4) ires=ires+2
3977 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3980 dc(j,ires)=sccor(j,iii)
3983 call sccenter(ires,iii,sccor)
3987 else if (card(:3).eq.'BRA') then
3991 itype(ires,molecule)=ntyp1_molec(molecule)-1
3992 nres_molec(molecule)=nres_molec(molecule)+1
3996 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3997 ! Fish out the ATOM cards.
3998 ! write(iout,*) 'card',card(1:20)
3999 ! print *,"ATU ",card(1:6), CARD(3:6)
4001 if (index(card(1:4),'ATOM').gt.0) then
4002 read (card(12:16),*) atom
4003 ! write (iout,*) "! ",atom," !",ires
4004 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
4005 if (card(14:16).eq.'LIP') then
4010 nres_molec(molecule)=nres_molec(molecule)+1
4011 itype(ires,molecule)=ntyp1_molec(molecule)
4020 nres_molec(molecule)=nres_molec(molecule)+1
4021 read (card(18:20),'(a3)') res
4023 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4024 if (UNRES_PDB) then!
4025 itype(ires,molecule)=rescode(ires,res,0,molecule)
4027 itype(ires,molecule)=rescode_lip(res,ires)
4030 read (card(23:26),*) ires
4031 read (card(18:20),'(a3)') res
4032 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
4033 ! & " ires_old",ires_old
4034 ! write (iout,*) "ishift",ishift," ishift1",ishift1
4035 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
4036 if (ires-ishift+ishift1.ne.ires_old) then
4037 ! Calculate the CM of the preceding residue.
4038 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
4040 ! write (iout,*) "Calculating sidechain center iii",iii
4043 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
4046 call sccenter(ires_old,iii,sccor)
4050 ! Start new residue.
4051 if (res.eq.'Cl-' .or. res.eq.'Na+') then
4054 else if (ibeg.eq.1) then
4055 write (iout,*) "BEG ires",ires
4057 if (res.ne.'GLY' .and. res.ne. 'ACE') then
4060 nres_molec(molecule)=nres_molec(molecule)+1
4062 ires=ires-ishift+ishift1
4064 ! write (iout,*) "ishift",ishift," ires",ires,&
4065 ! " ires_old",ires_old
4067 else if (ibeg.eq.2) then
4069 ishift=-ires_old+ires-1 !!!!!
4070 ishift1=ishift1-1 !!!!!
4071 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
4072 ires=ires-ishift+ishift1
4073 ! print *,ires,ishift,ishift1
4077 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
4078 ires=ires-ishift+ishift1
4081 ! print *,'atom',ires,atom
4082 if (res.eq.'ACE' .or. res.eq.'NHE') then
4085 if (atom.eq.'CA '.or.atom.eq.'N ') then
4087 itype(ires,molecule)=rescode(ires,res,0,molecule)
4089 ! nres_molec(molecule)=nres_molec(molecule)+1
4093 itype(ires,molecule)=rescode(ires,res2,0,molecule)
4094 ! nres_molec(molecule)=nres_molec(molecule)+1
4095 read (card(19:19),'(a1)') sugar
4096 isugar=sugarcode(sugar,ires)
4097 ! if (ibeg.eq.1) then
4101 ! print *,"ires=",ires,istype(ires)
4107 ires=ires-ishift+ishift1
4109 ! write (iout,*) "ires_old",ires_old," ires",ires
4110 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
4113 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
4114 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
4115 res.eq.'NHE'.and.atom(:2).eq.'HN') then
4116 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4117 ! print *,ires,ishift,ishift1
4118 ! write (iout,*) "backbone ",atom
4120 write (iout,'(2i3,2x,a,3f8.3)') &
4121 ires,itype(ires,1),res,(c(j,ires),j=1,3)
4124 nres_molec(molecule)=nres_molec(molecule)+1
4126 sccor(j,iii)=c(j,ires)
4128 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
4129 atom.eq."C2'" .or. atom.eq."C3'" &
4130 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
4131 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
4132 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
4133 ! print *,ires,ishift,ishift1
4137 ! sccor(j,iii)=c(j,ires)
4140 c(j,ires)=c(j,ires)+ccc(j)/5.0
4142 print *,counter,molecule
4143 if (counter.eq.5) then
4145 nres_molec(molecule)=nres_molec(molecule)+1
4148 ! sccor(j,iii)=c(j,ires)
4152 ! print *, "ATOM",atom(1:3)
4153 else if (atom.eq."C5'") then
4154 read (card(19:19),'(a1)') sugar
4155 isugar=sugarcode(sugar,ires)
4160 ! print *,ires,istype(ires)
4164 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
4165 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4166 nres_molec(molecule)=nres_molec(molecule)+1
4167 print *,"nres_molec(molecule)",nres_molec(molecule),ires
4171 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4173 else if ((atom.eq."C1'").and.unres_pdb) then
4175 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4176 ! write (*,*) card(23:27),ires,itype(ires,1)
4177 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
4178 atom.ne.'N' .and. atom.ne.'C' .and. &
4179 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
4180 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
4181 .and. atom.ne.'P '.and. &
4182 atom(1:1).ne.'H' .and. &
4183 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
4185 ! write (iout,*) "sidechain ",atom
4186 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4187 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
4188 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
4190 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
4194 ! print *,"IONS",ions,card(1:6)
4195 else if ((ions).and.(card(1:6).eq.'HETATM')) then
4196 if (firstion.eq.0) then
4200 dc(j,ires)=sccor(j,iii)
4203 call sccenter(ires,iii,sccor)
4206 read (card(12:16),*) atom
4207 print *,"HETATOM", atom(1:2)
4208 read (card(18:20),'(a3)') res
4209 if (atom(3:3).eq.'H') cycle
4210 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
4211 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
4212 .or.(atom(1:2).eq.'K ').or.(atom(1:2).eq.'ZN').or.&
4213 (atom(1:2).eq.'O ')) then
4215 print *,"I have water"
4216 if (molecule.ne.5) molecprev=molecule
4218 nres_molec(molecule)=nres_molec(molecule)+1
4219 print *,"HERE",nres_molec(molecule)
4220 if (res.ne.'WAT') res=res(2:3)//' '
4221 itype(ires,molecule)=rescode(ires,res,0,molecule)
4222 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
4226 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
4227 if (ires.eq.0) return
4228 ! Calculate dummy residue coordinates inside the "chain" of a multichain
4231 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
4233 nres_molec(molecule)=nres_molec(molecule)-2
4234 print *,'I have',nres, nres_molec(:)
4236 do k=1,4 ! ions are without dummy
4237 if (nres_molec(k).eq.0) cycle
4239 ! write (iout,*) i,itype(i,1)
4240 ! if (itype(i,1).eq.ntyp1) then
4241 ! write (iout,*) "dummy",i,itype(i,1)
4243 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
4244 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
4248 if (itype(i,k).eq.ntyp1_molec(k)) then
4249 if (itype(i+1,k).eq.ntyp1_molec(k)) then
4250 if (itype(i+2,k).eq.0) then
4251 ! print *,"masz sieczke"
4253 if (itype(i+2,j).ne.0) then
4255 itype(i+1,j)=ntyp1_molec(j)
4256 nres_molec(k)=nres_molec(k)-1
4257 nres_molec(j)=nres_molec(j)+1
4263 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
4264 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
4265 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
4266 ! if (unres_pdb) then
4267 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4268 ! print *,i,'tu dochodze'
4269 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
4277 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
4281 dcj=(c(j,i-2)-c(j,i-3))/2.0
4282 if (dcj.eq.0) dcj=1.23591524223
4287 else !itype(i+1,1).eq.ntyp1
4288 ! if (unres_pdb) then
4289 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4290 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
4297 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
4298 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
4302 dcj=(c(j,i+3)-c(j,i+2))/2.0
4303 if (dcj.eq.0) dcj=1.23591524223
4308 endif !itype(i+1,1).eq.ntyp1
4309 endif !itype.eq.ntyp1
4313 ! Calculate the CM of the last side chain.
4317 dc(j,ires)=sccor(j,iii)
4320 call sccenter(ires,iii,sccor)
4326 ! print *,"molecule",molecule
4327 if ((itype(nres,1).ne.10)) then
4330 if (molecule.eq.5) molecule=molecprev
4331 itype(nres,molecule)=ntyp1_molec(molecule)
4332 nres_molec(molecule)=nres_molec(molecule)+1
4333 ! if (unres_pdb) then
4334 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
4335 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
4342 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
4346 dcj=(c(j,nres-nres_molec(5)-2)-c(j,nres-nres_molec(5)-3))/2.0
4347 c(j,nres)=c(j,nres-nres_molec(5)-1)+dcj
4348 c(j,2*nres)=c(j,nres)
4352 ! print *,'I have',nres, nres_molec(:)
4354 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
4356 if (nres.ne.nres0) then
4357 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
4359 stop "Error nres value in WHAM input"
4362 !---------------------------------
4363 !el reallocate tables
4366 ! hfrag_alloc(j,i)=hfrag(j,i)
4369 ! bfrag_alloc(j,i)=bfrag(j,i)
4375 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
4376 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
4377 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
4378 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
4382 ! hfrag(j,i)=hfrag_alloc(j,i)
4387 ! bfrag(j,i)=bfrag_alloc(j,i)
4390 !el end reallocate tables
4391 !---------------------------------
4399 c(j,2*nres)=c(j,nres)
4402 if (itype(1,1).eq.ntyp1) then
4406 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4407 call refsys(2,3,4,e1,e2,e3,fail)
4414 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4415 c(j,1)=c(j,2)+1.9d0*(e1(j)-e2(j))/sqrt(2.0d0)
4419 dcj=(c(j,4)-c(j,3))/2.0
4425 ! First lets assign correct dummy to correct type of chain
4427 if (itype(1,1).eq.ntyp1) then
4428 if (itype(2,1).eq.0) then
4430 if (itype(2,j).ne.0) then
4432 itype(1,j)=ntyp1_molec(j)
4433 nres_molec(1)=nres_molec(1)-1
4434 nres_molec(j)=nres_molec(j)+1
4441 print *,'I have',nres, nres_molec(:)
4443 ! Copy the coordinates to reference coordinates
4449 ! Calculate internal coordinates.
4451 write (iout,'(/a)') &
4452 "Cartesian coordinates of the reference structure"
4453 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4454 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4456 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4457 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4458 (c(j,ires+nres),j=1,3)
4461 ! znamy już nres wiec mozna alokowac tablice
4462 ! Calculate internal coordinates.
4463 if(me.eq.king.or..not.out1file)then
4464 write (iout,'(a)') &
4465 "Backbone and SC coordinates as read from the PDB"
4467 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
4468 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
4469 (c(j,nres+ires),j=1,3)
4472 ! NOW LETS ROCK! SORTING
4473 allocate(c_temporary(3,2*nres))
4474 allocate(itype_temporary(nres,5))
4475 if (.not.allocated(molnum)) allocate(molnum(nres+1))
4476 if (.not.allocated(istype)) write(iout,*) &
4477 "SOMETHING WRONG WITH ISTYTPE"
4478 allocate(istype_temp(nres))
4479 itype_temporary(:,:)=0
4483 if (itype(i,k).ne.0) then
4485 c_temporary(j,seqalingbegin)=c(j,i)
4486 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
4488 itype_temporary(seqalingbegin,k)=itype(i,k)
4489 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
4490 istype_temp(seqalingbegin)=istype(i)
4491 molnum(seqalingbegin)=k
4492 seqalingbegin=seqalingbegin+1
4498 c(j,i)=c_temporary(j,i)
4500 if ((molnum(i-nres).eq.4)) c(j,i)=0.0d0
4506 itype(i,k)=itype_temporary(i,k)
4507 istype(i)=istype_temp(i)
4510 if ((itype(1,1).eq.ntyp1).and.itype(2,5).ne.0) then
4511 ! I have only ions now dummy atoms in the system
4513 itype(1,5)=itype(2,5)
4519 itype(i,5)=itype(i+1,5)
4526 nres_molec(1)=nres_molec(1)-1
4528 ! if (itype(1,1).eq.ntyp1) then
4531 ! if (unres_pdb) then
4532 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
4533 ! call refsys(2,3,4,e1,e2,e3,fail)
4540 ! c(j,1)=c(j,2)-1.9d0*e2(j)
4544 ! dcj=(c(j,4)-c(j,3))/2.0
4546 ! c(j,nres+1)=c(j,1)
4552 write (iout,'(/a)') &
4553 "Cartesian coordinates of the reference structure after sorting"
4554 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
4555 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
4557 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
4558 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
4559 (c(j,ires+nres),j=1,3)
4563 print *,seqalingbegin,nres
4564 if(.not.allocated(vbld)) then
4565 allocate(vbld(2*nres))
4570 if(.not.allocated(vbld_inv)) then
4571 allocate(vbld_inv(2*nres))
4577 if(.not.allocated(theta)) then
4578 allocate(theta(nres+2))
4582 if(.not.allocated(phi)) allocate(phi(nres+2))
4583 if(.not.allocated(alph)) allocate(alph(nres+2))
4584 if(.not.allocated(omeg)) allocate(omeg(nres+2))
4585 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
4586 if(.not.allocated(phiref)) allocate(phiref(nres+2))
4587 if(.not.allocated(costtab)) allocate(costtab(nres))
4588 if(.not.allocated(sinttab)) allocate(sinttab(nres))
4589 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
4590 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
4591 if(.not.allocated(xxref)) allocate(xxref(nres))
4592 if(.not.allocated(yyref)) allocate(yyref(nres))
4593 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
4594 if(.not.allocated(dc_norm)) then
4595 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
4596 allocate(dc_norm(3,0:2*nres+2))
4599 write(iout,*) "before int_from_cart"
4600 call int_from_cart(.true.,.false.)
4601 call sc_loc_geom(.false.)
4602 write(iout,*) "after int_from_cart"
4606 thetaref(i)=theta(i)
4609 write(iout,*) "after thetaref"
4617 dc(j,i)=c(j,i+1)-c(j,i)
4618 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
4623 dc(j,i+nres)=c(j,i+nres)-c(j,i)
4624 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
4626 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
4630 ! Copy the coordinates to reference coordinates
4631 ! Splits to single chain if occurs
4635 ! cref(j,i,cou)=c(j,i)
4641 ! chomo(j,i,k)=c(j,i)
4644 ! write(iout,*) "after chomo"
4646 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4647 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4648 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4649 !-----------------------------
4653 write (iout,*) "symetr", symetr
4656 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4658 ! if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4659 ! chain_length=lll-1
4661 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4666 cref(j,i,cou)=c(j,i)
4667 cref(j,i+nres,cou)=c(j,i+nres)
4669 chain_rep(j,lll,kkk)=c(j,i)
4670 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4675 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4677 ! do kkk=1,chain_length
4678 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4682 ! makes copy of chains
4683 write (iout,*) "symetr", symetr
4688 ! if (symetr.gt.1) then
4689 ! call permut(symetr)
4695 ! write(iout,*) (tabperm(i,kkk),kkk=1,4)
4700 ! icha=tabperm(i,kkk)
4701 ! write (iout,*) i,icha
4702 ! do lll=1,chain_length
4704 ! if (cou.le.nres) then
4706 ! kupa=mod(lll,chain_length)
4707 ! iprzes=(kkk-1)*chain_length+lll
4708 ! if (kupa.eq.0) kupa=chain_length
4709 ! write (iout,*) "kupa", kupa
4710 ! cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4711 ! cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4718 !-koniec robienia kopii
4721 write (iout,*) "nowa struktura", nperm
4723 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4725 cref(3,i,kkk),cref(1,nres+i,kkk),&
4726 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4728 100 format (//' alpha-carbon coordinates ',&
4729 ' centroid coordinates'/ &
4730 ' ', 6X,'X',11X,'Y',11X,'Z', &
4731 10X,'X',11X,'Y',11X,'Z')
4732 110 format (a,'(',i5,')',6f12.5)
4738 bfrag(i,j)=bfrag(i,j)-ishift
4744 hfrag(i,j)=hfrag(i,j)-ishift
4750 end subroutine readpdb
4751 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4752 !-----------------------------------------------------------------------------
4754 !-----------------------------------------------------------------------------
4755 subroutine read_control
4769 use random, only: random_init
4770 ! implicit real*8 (a-h,o-z)
4771 ! include 'DIMENSIONS'
4773 use prng, only:prng_restart
4775 logical :: OKRandom!, prng_restart
4778 ! include 'COMMON.IOUNITS'
4779 ! include 'COMMON.TIME1'
4780 ! include 'COMMON.THREAD'
4781 ! include 'COMMON.SBRIDGE'
4782 ! include 'COMMON.CONTROL'
4783 ! include 'COMMON.MCM'
4784 ! include 'COMMON.MAP'
4785 ! include 'COMMON.HEADER'
4786 ! include 'COMMON.CSA'
4787 ! include 'COMMON.CHAIN'
4788 ! include 'COMMON.MUCA'
4789 ! include 'COMMON.MD'
4790 ! include 'COMMON.FFIELD'
4791 ! include 'COMMON.INTERACT'
4792 ! include 'COMMON.SETUP'
4793 !el integer :: KDIAG,ICORFL,IXDR
4794 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4795 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4796 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4797 ! character(len=80) :: ucase
4798 character(len=640) :: controlcard
4800 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4802 usampl=.false. ! the default value of usample should be 0
4803 ! write(iout,*) "KURWA2",usampl
4807 read (INP,'(a)') titel
4808 call card_concat(controlcard,.true.)
4809 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4810 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4811 call reada(controlcard,'SEED',seed,0.0D0)
4812 call random_init(seed)
4813 ! Set up the time limit (caution! The time must be input in minutes!)
4814 read_cart=index(controlcard,'READ_CART').gt.0
4815 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4816 call readi(controlcard,'SYM',symetr,1)
4817 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4818 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4819 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4820 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4821 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4822 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4823 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4824 call reada(controlcard,'DRMS',drms,0.1D0)
4825 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
4826 read_homol_frag = index(controlcard,"READ_HOMOL_FRAG").gt.0
4827 out_template_coord = index(controlcard,"OUT_TEMPLATE_COORD").gt.0
4828 out_template_restr = index(controlcard,"OUT_TEMPLATE_RESTR").gt.0
4829 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4830 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4831 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4832 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4833 write (iout,'(a,f10.1)')'DRMS = ',drms
4834 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4835 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4837 call readi(controlcard,'NZ_START',nz_start,0)
4838 call readi(controlcard,'NZ_END',nz_end,0)
4839 call readi(controlcard,'IZ_SC',iz_sc,0)
4840 timlim=60.0D0*timlim
4841 safety = 60.0d0*safety
4844 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4845 !C SHIELD keyword sets if the shielding effect of side-chains is used
4846 !C 0 denots no shielding is used all peptide are equally despite the
4847 !C solvent accesible area
4848 !C 1 the newly introduced function
4849 !C 2 reseved for further possible developement
4850 call readi(controlcard,'SHIELD',shield_mode,0)
4851 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4852 write(iout,*) "shield_mode",shield_mode
4853 call readi(controlcard,'TORMODE',tor_mode,0)
4854 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4855 write(iout,*) "torsional and valence angle mode",tor_mode
4857 !C Varibles set size of box
4858 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4859 protein=index(controlcard,"PROTEIN").gt.0
4860 ions=index(controlcard,"IONS").gt.0
4861 fodson=index(controlcard,"FODSON").gt.0
4862 call readi(controlcard,'OLDION',oldion,1)
4863 nucleic=index(controlcard,"NUCLEIC").gt.0
4864 write (iout,*) "with_theta_constr ",with_theta_constr
4865 AFMlog=(index(controlcard,'AFM'))
4866 selfguide=(index(controlcard,'SELFGUIDE'))
4867 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4868 call readi(controlcard,'GENCONSTR',genconstr,0)
4869 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4870 call reada(controlcard,'BOXY',boxysize,100.0d0)
4871 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4872 call readi(controlcard,'TUBEMOD',tubemode,0)
4873 print *,"SCELE",scelemode
4874 call readi(controlcard,"SCELEMODE",scelemode,0)
4875 print *,"SCELE",scelemode
4877 ! elemode = 0 is orignal UNRES electrostatics
4878 ! elemode = 1 is "Momo" potentials in progress
4879 ! elemode = 2 is in development EVALD
4882 write (iout,*) TUBEmode,"TUBEMODE"
4883 if (TUBEmode.gt.0) then
4884 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4885 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4886 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4887 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4888 call reada(controlcard,"VNANO",velnanoconst,0.0d0)
4889 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4890 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4891 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4892 buftubebot=bordtubebot+tubebufthick
4893 buftubetop=bordtubetop-tubebufthick
4896 ! CUTOFFF ON ELECTROSTATICS
4897 call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
4898 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4900 write(iout,*) "R_CUT_ELE=",r_cut_ele,rlamb_ele
4901 call reada(controlcard,"R_CUT_MART",r_cut_mart,15.0d0)
4902 call reada(controlcard,"LAMBDA_MART",rlamb_mart,0.3d0)
4903 call reada(controlcard,"R_CUT_ANG",r_cut_ang,4.2d0)
4905 call reada(controlcard,"DELTA",graddelta,1.0d-5)
4906 ! Lipidec parameters
4907 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4908 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4909 if (lipthick.gt.0.0d0) then
4910 bordliptop=(boxzsize+lipthick)/2.0
4911 bordlipbot=bordliptop-lipthick
4912 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4913 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4914 buflipbot=bordlipbot+lipbufthick
4915 bufliptop=bordliptop-lipbufthick
4916 if ((lipbufthick*2.0d0).gt.lipthick) &
4917 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4918 endif !lipthick.gt.0
4919 write(iout,*) "bordliptop=",bordliptop
4920 write(iout,*) "bordlipbot=",bordlipbot
4921 write(iout,*) "bufliptop=",bufliptop
4922 write(iout,*) "buflipbot=",buflipbot
4923 write (iout,*) "SHIELD MODE",shield_mode
4925 !C-------------------------
4926 minim=(index(controlcard,'MINIMIZE').gt.0)
4927 dccart=(index(controlcard,'CART').gt.0)
4928 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4929 overlapsc=.not.overlapsc
4930 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4931 searchsc=.not.searchsc
4932 sideadd=(index(controlcard,'SIDEADD').gt.0)
4933 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4934 outpdb=(index(controlcard,'PDBOUT').gt.0)
4935 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4936 pdbref=(index(controlcard,'PDBREF').gt.0)
4937 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4938 indpdb=index(controlcard,'PDBSTART')
4939 extconf=(index(controlcard,'EXTCONF').gt.0)
4940 call readi(controlcard,'IPRINT',iprint,0)
4941 call readi(controlcard,'MAXGEN',maxgen,10000)
4942 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4943 call readi(controlcard,"KDIAG",kdiag,0)
4944 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4945 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4946 write (iout,*) "RESCALE_MODE",rescale_mode
4947 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4948 if (index(controlcard,'REGULAR').gt.0.0D0) then
4949 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4953 if (index(controlcard,'CHECKGRAD').gt.0) then
4955 if (index(controlcard,'CART').gt.0) then
4957 elseif (index(controlcard,'CARINT').gt.0) then
4962 elseif (index(controlcard,'THREAD').gt.0) then
4964 call readi(controlcard,'THREAD',nthread,0)
4965 if (nthread.gt.0) then
4966 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4969 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4970 stop 'Error termination in Read_Control.'
4972 else if (index(controlcard,'MCMA').gt.0) then
4974 else if (index(controlcard,'MCEE').gt.0) then
4976 else if (index(controlcard,'MULTCONF').gt.0) then
4978 else if (index(controlcard,'MAP').gt.0) then
4980 call readi(controlcard,'MAP',nmap,0)
4981 else if (index(controlcard,'CSA').gt.0) then
4983 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4985 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4988 !fcm else if (index(controlcard,'MCMF').gt.0) then
4990 else if (index(controlcard,'SOFTREG').gt.0) then
4992 else if (index(controlcard,'CHECK_BOND').gt.0) then
4994 else if (index(controlcard,'TEST').gt.0) then
4996 else if (index(controlcard,'MD').gt.0) then
4998 else if (index(controlcard,'RE ').gt.0) then
5002 lmuca=index(controlcard,'MUCA').gt.0
5003 call readi(controlcard,'MUCADYN',mucadyn,0)
5004 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
5005 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
5007 write (iout,*) 'MUCADYN=',mucadyn
5008 write (iout,*) 'MUCASMOOTH=',muca_smooth
5011 iscode=index(controlcard,'ONE_LETTER')
5012 indphi=index(controlcard,'PHI')
5013 indback=index(controlcard,'BACK')
5014 iranconf=index(controlcard,'RAND_CONF')
5015 i2ndstr=index(controlcard,'USE_SEC_PRED')
5016 gradout=index(controlcard,'GRADOUT').gt.0
5017 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
5018 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
5019 if (me.eq.king .or. .not.out1file ) &
5020 write (iout,*) "DISTCHAINMAX",distchainmax
5022 if(me.eq.king.or..not.out1file) &
5023 write (iout,'(2a)') diagmeth(kdiag),&
5024 ' routine used to diagonalize matrices.'
5025 if (shield_mode.gt.0) then
5026 pi=4.0D0*datan(1.0D0)
5027 !C VSolvSphere the volume of solving sphere
5029 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
5030 !C there will be no distinction between proline peptide group and normal peptide
5031 !C group in case of shielding parameters
5032 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
5033 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
5034 write (iout,*) VSolvSphere,VSolvSphere_div
5035 !C long axis of side chain
5037 ! long_r_sidechain(i)=vbldsc0(1,i)
5038 ! short_r_sidechain(i)=sigma0(i)
5039 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
5046 end subroutine read_control
5047 !-----------------------------------------------------------------------------
5048 subroutine read_REMDpar
5050 ! Read REMD settings
5057 use control_data, only:out1file
5058 ! implicit real*8 (a-h,o-z)
5059 ! include 'DIMENSIONS'
5060 ! include 'COMMON.IOUNITS'
5061 ! include 'COMMON.TIME1'
5062 ! include 'COMMON.MD'
5065 !el include 'COMMON.LANGEVIN'
5067 !el include 'COMMON.LANGEVIN.lang0'
5069 ! include 'COMMON.INTERACT'
5070 ! include 'COMMON.NAMES'
5071 ! include 'COMMON.GEO'
5072 ! include 'COMMON.REMD'
5073 ! include 'COMMON.CONTROL'
5074 ! include 'COMMON.SETUP'
5075 ! character(len=80) :: ucase
5076 character(len=320) :: controlcard
5077 character(len=3200) :: controlcard1
5078 integer :: iremd_m_total
5081 ! real(kind=8) :: var,ene
5083 if(me.eq.king.or..not.out1file) &
5084 write (iout,*) "REMD setup"
5086 call card_concat(controlcard,.true.)
5087 call readi(controlcard,"NREP",nrep,3)
5088 call readi(controlcard,"NSTEX",nstex,1000)
5089 call reada(controlcard,"RETMIN",retmin,10.0d0)
5090 call reada(controlcard,"RETMAX",retmax,1000.0d0)
5091 mremdsync=(index(controlcard,'SYNC').gt.0)
5092 call readi(controlcard,"NSYN",i_sync_step,100)
5093 restart1file=(index(controlcard,'REST1FILE').gt.0)
5094 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
5095 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
5096 if(max_cache_traj_use.gt.max_cache_traj) &
5097 max_cache_traj_use=max_cache_traj
5098 if(me.eq.king.or..not.out1file) then
5099 !d if (traj1file) then
5100 !rc caching is in testing - NTWX is not ignored
5101 !d write (iout,*) "NTWX value is ignored"
5102 !d write (iout,*) " trajectory is stored to one file by master"
5103 !d write (iout,*) " before exchange at NSTEX intervals"
5105 write (iout,*) "NREP= ",nrep
5106 write (iout,*) "NSTEX= ",nstex
5107 write (iout,*) "SYNC= ",mremdsync
5108 write (iout,*) "NSYN= ",i_sync_step
5109 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
5112 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
5113 if (index(controlcard,'TLIST').gt.0) then
5115 call card_concat(controlcard1,.true.)
5116 read(controlcard1,*) (remd_t(i),i=1,nrep)
5117 if(me.eq.king.or..not.out1file) &
5118 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
5121 if (index(controlcard,'MLIST').gt.0) then
5123 call card_concat(controlcard1,.true.)
5124 read(controlcard1,*) (remd_m(i),i=1,nrep)
5125 if(me.eq.king.or..not.out1file) then
5126 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
5129 iremd_m_total=iremd_m_total+remd_m(i)
5131 write (iout,*) 'Total number of replicas ',iremd_m_total
5134 if(me.eq.king.or..not.out1file) &
5135 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
5137 end subroutine read_REMDpar
5138 !-----------------------------------------------------------------------------
5139 subroutine read_MDpar
5143 use control_data, only: r_cut,rlamb,out1file,r_cut_mart,rlamb_mart
5145 use geometry_data, only: pi
5147 ! implicit real*8 (a-h,o-z)
5148 ! include 'DIMENSIONS'
5149 ! include 'COMMON.IOUNITS'
5150 ! include 'COMMON.TIME1'
5151 ! include 'COMMON.MD'
5154 !el include 'COMMON.LANGEVIN'
5156 !el include 'COMMON.LANGEVIN.lang0'
5158 ! include 'COMMON.INTERACT'
5159 ! include 'COMMON.NAMES'
5160 ! include 'COMMON.GEO'
5161 ! include 'COMMON.SETUP'
5162 ! include 'COMMON.CONTROL'
5163 ! include 'COMMON.SPLITELE'
5164 ! character(len=80) :: ucase
5165 character(len=320) :: controlcard
5170 call card_concat(controlcard,.true.)
5171 call readi(controlcard,"NSTEP",n_timestep,1000000)
5172 call readi(controlcard,"NTWE",ntwe,100)
5173 call readi(controlcard,"NTWX",ntwx,1000)
5174 call readi(controlcard,"NFOD",nfodstep,100)
5175 call reada(controlcard,"DT",d_time,1.0d-1)
5176 call reada(controlcard,"DVMAX",dvmax,2.0d1)
5177 call reada(controlcard,"DAMAX",damax,1.0d1)
5178 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
5179 call readi(controlcard,"LANG",lang,0)
5180 RESPA = index(controlcard,"RESPA") .gt. 0
5181 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
5182 ntime_split0=ntime_split
5183 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
5184 ntime_split0=ntime_split
5185 call reada(controlcard,"R_CUT",r_cut,2.0d0)
5186 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
5187 rest = index(controlcard,"REST").gt.0
5188 tbf = index(controlcard,"TBF").gt.0
5189 usampl = index(controlcard,"USAMPL").gt.0
5190 ! write(iout,*) "KURWA",usampl
5191 mdpdb = index(controlcard,"MDPDB").gt.0
5192 call reada(controlcard,"T_BATH",t_bath,300.0d0)
5193 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
5194 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
5195 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
5196 if (count_reset_moment.eq.0) count_reset_moment=1000000000
5197 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
5198 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
5199 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
5200 if (count_reset_vel.eq.0) count_reset_vel=1000000000
5201 large = index(controlcard,"LARGE").gt.0
5202 print_compon = index(controlcard,"PRINT_COMPON").gt.0
5203 rattle = index(controlcard,"RATTLE").gt.0
5204 preminim=(index(controlcard,'PREMINIM').gt.0)
5205 forceminim=(index(controlcard,'FORCEMINIM').gt.0)
5206 write (iout,*) "PREMINIM ",preminim
5207 dccart=(index(controlcard,'CART').gt.0)
5208 if (preminim) call read_minim
5209 ! if performing umbrella sampling, fragments constrained are read from the fragment file
5215 if(me.eq.king.or..not.out1file) then
5217 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
5219 write (iout,'(a)') "The units are:"
5220 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
5221 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
5222 " acceleration: angstrom/(48.9 fs)**2"
5223 write (iout,'(a)') "energy: kcal/mol, temperature: K"
5225 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
5226 write (iout,'(a60,f10.5,a)') &
5227 "Initial time step of numerical integration:",d_time,&
5229 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
5231 write (iout,'(2a,i4,a)') &
5232 "A-MTS algorithm used; initial time step for fast-varying",&
5233 " short-range forces split into",ntime_split," steps."
5234 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
5235 r_cut," lambda",rlamb
5237 write (iout,'(2a,f10.5)') &
5238 "Maximum acceleration threshold to reduce the time step",&
5239 "/increase split number:",damax
5240 write (iout,'(2a,f10.5)') &
5241 "Maximum predicted energy drift to reduce the timestep",&
5242 "/increase split number:",edriftmax
5243 write (iout,'(a60,f10.5)') &
5244 "Maximum velocity threshold to reduce velocities:",dvmax
5245 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
5246 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
5247 if (rattle) write (iout,'(a60)') &
5248 "Rattle algorithm used to constrain the virtual bonds"
5252 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
5253 call reada(controlcard,"RWAT",rwat,1.4d0)
5254 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
5255 surfarea=index(controlcard,"SURFAREA").gt.0
5256 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
5257 if(me.eq.king.or..not.out1file)then
5258 write (iout,'(/a,$)') "Langevin dynamics calculation"
5260 write (iout,'(a/)') &
5261 " with direct integration of Langevin equations"
5262 else if (lang.eq.2) then
5263 write (iout,'(a/)') " with TINKER stochasic MD integrator"
5264 else if (lang.eq.3) then
5265 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
5266 else if (lang.eq.4) then
5267 write (iout,'(a/)') " in overdamped mode"
5269 write (iout,'(//a,i5)') &
5270 "=========== ERROR: Unknown Langevin dynamics mode:",lang
5273 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5274 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
5275 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
5276 write (iout,'(a60,f10.5)') &
5277 "Scaling factor of the friction forces:",scal_fric
5278 if (surfarea) write (iout,'(2a,i10,a)') &
5279 "Friction coefficients will be scaled by solvent-accessible",&
5280 " surface area every",reset_fricmat," steps."
5282 ! Calculate friction coefficients and bounds of stochastic forces
5283 eta=6*pi*cPoise*etawat
5284 if(me.eq.king.or..not.out1file) &
5285 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
5288 do j=1,5 !types of molecules
5289 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
5290 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
5292 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
5293 do j=1,5 !types of molecules
5295 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
5296 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
5300 if(me.eq.king.or..not.out1file)then
5302 write (iout,'(/2a/)') &
5303 "Radii of site types and friction coefficients and std's of",&
5304 " stochastic forces of fully exposed sites"
5305 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(j),stdfp*dsqrt(gamp(j))
5308 write (iout,'(a5,f5.2,2f10.5)') restyp(i,j),restok(i,j),&
5309 gamsc(i,j),stdfsc(i,j)*dsqrt(gamsc(i,j))
5314 if(me.eq.king.or..not.out1file)then
5315 write (iout,'(a)') "Berendsen bath calculation"
5316 write (iout,'(a60,f10.5)') "Temperature:",t_bath
5317 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
5319 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
5320 count_reset_moment," steps"
5322 write (iout,'(a,i10,a)') &
5323 "Velocities will be reset at random every",count_reset_vel,&
5327 if(me.eq.king.or..not.out1file) &
5328 write (iout,'(a31)') "Microcanonical mode calculation"
5330 if(me.eq.king.or..not.out1file)then
5331 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
5333 write(iout,*) "MD running with constraints."
5334 write(iout,*) "Equilibration time ", eq_time, " mtus."
5335 write(iout,*) "Constraining ", nfrag," fragments."
5336 write(iout,*) "Length of each fragment, weight and q0:"
5338 write (iout,*) "Set of restraints #",iset
5340 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
5341 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
5343 write(iout,*) "constraints between ", npair, "fragments."
5344 write(iout,*) "constraint pairs, weights and q0:"
5346 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
5347 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
5349 write(iout,*) "angle constraints within ", nfrag_back,&
5350 "backbone fragments."
5351 write(iout,*) "fragment, weights:"
5353 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
5354 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
5355 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
5358 iset=mod(kolor,nset)+1
5361 if(me.eq.king.or..not.out1file) &
5362 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
5364 end subroutine read_MDpar
5365 !-----------------------------------------------------------------------------
5369 ! implicit real*8 (a-h,o-z)
5370 ! include 'DIMENSIONS'
5371 ! include 'COMMON.MAP'
5372 ! include 'COMMON.IOUNITS'
5373 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
5374 character(len=80) :: mapcard !,ucase
5377 ! real(kind=8) :: var,ene
5380 read (inp,'(a)') mapcard
5381 mapcard=ucase(mapcard)
5382 if (index(mapcard,'PHI').gt.0) then
5384 else if (index(mapcard,'THE').gt.0) then
5386 else if (index(mapcard,'ALP').gt.0) then
5388 else if (index(mapcard,'OME').gt.0) then
5391 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
5392 stop 'Error - illegal variable spec in MAP card.'
5394 call readi (mapcard,'RES1',res1(imap),0)
5395 call readi (mapcard,'RES2',res2(imap),0)
5396 if (res1(imap).eq.0) then
5397 res1(imap)=res2(imap)
5398 else if (res2(imap).eq.0) then
5399 res2(imap)=res1(imap)
5401 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
5402 write (iout,'(a)') &
5403 'Error - illegal definition of variable group in MAP.'
5404 stop 'Error - illegal definition of variable group in MAP.'
5406 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
5407 call reada(mapcard,'TO',ang_to(imap),0.0D0)
5408 call readi(mapcard,'NSTEP',nstep(imap),0)
5409 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
5410 write (iout,'(a)') &
5411 'Illegal boundary and/or step size specification in MAP.'
5412 stop 'Illegal boundary and/or step size specification in MAP.'
5416 end subroutine map_read
5417 !-----------------------------------------------------------------------------
5420 use control_data, only: vdisulf
5422 ! implicit real*8 (a-h,o-z)
5423 ! include 'DIMENSIONS'
5424 ! include 'COMMON.IOUNITS'
5425 ! include 'COMMON.GEO'
5426 ! include 'COMMON.CSA'
5427 ! include 'COMMON.BANK'
5428 ! include 'COMMON.CONTROL'
5429 ! character(len=80) :: ucase
5430 character(len=620) :: mcmcard
5432 ! integer :: ntf,ik,iw_pdb
5433 ! real(kind=8) :: var,ene
5435 call card_concat(mcmcard,.true.)
5437 call readi(mcmcard,'NCONF',nconf,50)
5438 call readi(mcmcard,'NADD',nadd,0)
5439 call readi(mcmcard,'JSTART',jstart,1)
5440 call readi(mcmcard,'JEND',jend,1)
5441 call readi(mcmcard,'NSTMAX',nstmax,500000)
5442 call readi(mcmcard,'N0',n0,1)
5443 call readi(mcmcard,'N1',n1,6)
5444 call readi(mcmcard,'N2',n2,4)
5445 call readi(mcmcard,'N3',n3,0)
5446 call readi(mcmcard,'N4',n4,0)
5447 call readi(mcmcard,'N5',n5,0)
5448 call readi(mcmcard,'N6',n6,10)
5449 call readi(mcmcard,'N7',n7,0)
5450 call readi(mcmcard,'N8',n8,0)
5451 call readi(mcmcard,'N9',n9,0)
5452 call readi(mcmcard,'N14',n14,0)
5453 call readi(mcmcard,'N15',n15,0)
5454 call readi(mcmcard,'N16',n16,0)
5455 call readi(mcmcard,'N17',n17,0)
5456 call readi(mcmcard,'N18',n18,0)
5458 vdisulf=(index(mcmcard,'DYNSS').gt.0)
5460 call readi(mcmcard,'NDIFF',ndiff,2)
5461 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
5462 call readi(mcmcard,'IS1',is1,1)
5463 call readi(mcmcard,'IS2',is2,8)
5464 call readi(mcmcard,'NRAN0',nran0,4)
5465 call readi(mcmcard,'NRAN1',nran1,2)
5466 call readi(mcmcard,'IRR',irr,1)
5467 call readi(mcmcard,'NSEED',nseed,20)
5468 call readi(mcmcard,'NTOTAL',ntotal,10000)
5469 call reada(mcmcard,'CUT1',cut1,2.0d0)
5470 call reada(mcmcard,'CUT2',cut2,5.0d0)
5471 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
5472 call readi(mcmcard,'ICMAX',icmax,3)
5473 call readi(mcmcard,'IRESTART',irestart,0)
5474 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
5477 call reada(mcmcard,'DELE',dele,20.0d0)
5478 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
5479 call readi(mcmcard,'IREF',iref,0)
5480 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
5481 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
5482 call readi(mcmcard,'NCONF_IN',nconf_in,0)
5483 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
5484 write (iout,*) "NCONF_IN",nconf_in
5486 end subroutine csaread
5487 !-----------------------------------------------------------------------------
5491 use control_data, only: MaxMoveType
5494 ! implicit real*8 (a-h,o-z)
5495 ! include 'DIMENSIONS'
5496 ! include 'COMMON.MCM'
5497 ! include 'COMMON.MCE'
5498 ! include 'COMMON.IOUNITS'
5499 ! character(len=80) :: ucase
5500 character(len=320) :: mcmcard
5503 ! real(kind=8) :: var,ene
5505 call card_concat(mcmcard,.true.)
5506 call readi(mcmcard,'MAXACC',maxacc,100)
5507 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
5508 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
5509 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
5510 call readi(mcmcard,'MAXREPM',maxrepm,200)
5511 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
5512 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
5513 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
5514 call reada(mcmcard,'E_UP',e_up,5.0D0)
5515 call reada(mcmcard,'DELTE',delte,0.1D0)
5516 call readi(mcmcard,'NSWEEP',nsweep,5)
5517 call readi(mcmcard,'NSTEPH',nsteph,0)
5518 call readi(mcmcard,'NSTEPC',nstepc,0)
5519 call reada(mcmcard,'TMIN',tmin,298.0D0)
5520 call reada(mcmcard,'TMAX',tmax,298.0D0)
5521 call readi(mcmcard,'NWINDOW',nwindow,0)
5522 call readi(mcmcard,'PRINT_MC',print_mc,0)
5523 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
5524 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
5525 ent_read=(index(mcmcard,'ENT_READ').gt.0)
5526 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
5527 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
5528 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
5529 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
5530 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
5531 if (nwindow.gt.0) then
5532 allocate(winstart(nwindow)) !!el (maxres)
5533 allocate(winend(nwindow)) !!el
5534 allocate(winlen(nwindow)) !!el
5535 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
5537 winlen(i)=winend(i)-winstart(i)+1
5540 if (tmax.lt.tmin) tmax=tmin
5541 if (tmax.eq.tmin) then
5545 if (nstepc.gt.0 .and. nsteph.gt.0) then
5546 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
5547 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
5549 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
5550 ! Probabilities of different move types
5551 sumpro_type(0)=0.0D0
5552 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
5553 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
5554 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
5555 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
5556 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
5557 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
5558 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
5560 print *,'i',i,' sumprotype',sumpro_type(i)
5561 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
5562 print *,'i',i,' sumprotype',sumpro_type(i)
5565 end subroutine mcmread
5566 !-----------------------------------------------------------------------------
5567 subroutine read_minim
5570 ! implicit real*8 (a-h,o-z)
5571 ! include 'DIMENSIONS'
5572 ! include 'COMMON.MINIM'
5573 ! include 'COMMON.IOUNITS'
5574 ! character(len=80) :: ucase
5575 character(len=320) :: minimcard
5577 ! integer :: ntf,ik,iw_pdb
5578 ! real(kind=8) :: var,ene
5580 call card_concat(minimcard,.true.)
5581 call readi(minimcard,'MAXMIN',maxmin,2000)
5582 call readi(minimcard,'MAXFUN',maxfun,5000)
5583 call readi(minimcard,'MINMIN',minmin,maxmin)
5584 call readi(minimcard,'MINFUN',minfun,maxmin)
5585 call reada(minimcard,'TOLF',tolf,1.0D-2)
5586 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
5587 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
5588 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
5589 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
5590 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
5591 'Options in energy minimization:'
5592 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
5593 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
5594 'MinMin:',MinMin,' MinFun:',MinFun,&
5595 ' TolF:',TolF,' RTolF:',RTolF
5597 end subroutine read_minim
5598 !-----------------------------------------------------------------------------
5599 subroutine openunits
5601 use MD_data, only: usampl
5604 use control_data, only:out1file
5605 use control, only: getenv_loc
5606 ! implicit real*8 (a-h,o-z)
5607 ! include 'DIMENSIONS'
5610 character(len=16) :: form,nodename
5611 integer :: nodelen,ierror,npos
5613 ! include 'COMMON.SETUP'
5614 ! include 'COMMON.IOUNITS'
5615 ! include 'COMMON.MD'
5616 ! include 'COMMON.CONTROL'
5617 integer :: lenpre,lenpot,lentmp !,ilen
5619 character(len=3) :: out1file_text !,ucase
5620 character(len=3) :: ll
5623 ! integer :: ntf,ik,iw_pdb
5624 ! real(kind=8) :: var,ene
5626 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
5627 call getenv_loc("PREFIX",prefix)
5629 call getenv_loc("POT",pot)
5630 call getenv_loc("DIRTMP",tmpdir)
5631 call getenv_loc("CURDIR",curdir)
5632 call getenv_loc("OUT1FILE",out1file_text)
5633 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
5634 out1file_text=ucase(out1file_text)
5635 if (out1file_text(1:1).eq."Y") then
5638 out1file=fg_rank.gt.0
5643 if (lentmp.gt.0) then
5644 write (*,'(80(1h!))')
5645 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
5646 write (*,'(80(1h!))')
5647 write (*,*)"All output files will be on node /tmp directory."
5649 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
5650 if (me.eq.king) then
5651 write (*,*) "The master node is ",nodename
5652 else if (fg_rank.eq.0) then
5653 write (*,*) "I am the CG slave node ",nodename
5655 write (*,*) "I am the FG slave node ",nodename
5658 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5659 lenpre = lentmp+lenpre+1
5661 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5662 ! Get the names and open the input files
5663 #if defined(WINIFL) || defined(WINPGI)
5664 open(1,file=pref_orig(:ilen(pref_orig))// &
5665 '.inp',status='old',readonly,shared)
5666 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5667 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5668 ! Get parameter filenames and open the parameter files.
5669 call getenv_loc('BONDPAR',bondname)
5670 open (ibond,file=bondname,status='old',readonly,shared)
5671 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5672 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5673 call getenv_loc('THETPAR',thetname)
5674 open (ithep,file=thetname,status='old',readonly,shared)
5675 call getenv_loc('ROTPAR',rotname)
5676 open (irotam,file=rotname,status='old',readonly,shared)
5677 call getenv_loc('TORPAR',torname)
5678 open (itorp,file=torname,status='old',readonly,shared)
5679 call getenv_loc('TORDPAR',tordname)
5680 open (itordp,file=tordname,status='old',readonly,shared)
5681 call getenv_loc('FOURIER',fouriername)
5682 open (ifourier,file=fouriername,status='old',readonly,shared)
5683 call getenv_loc('ELEPAR',elename)
5684 open (ielep,file=elename,status='old',readonly,shared)
5685 call getenv_loc('SIDEPAR',sidename)
5686 open (isidep,file=sidename,status='old',readonly,shared)
5688 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5689 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5690 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5691 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5692 call getenv_loc('TORPAR_NUCL',torname_nucl)
5693 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5694 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5695 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5696 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5697 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5700 #elif (defined CRAY) || (defined AIX)
5701 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5703 ! print *,"Processor",myrank," opened file 1"
5704 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5705 ! print *,"Processor",myrank," opened file 9"
5706 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5707 ! Get parameter filenames and open the parameter files.
5708 call getenv_loc('BONDPAR',bondname)
5709 open (ibond,file=bondname,status='old',action='read')
5710 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5711 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5713 ! print *,"Processor",myrank," opened file IBOND"
5714 call getenv_loc('THETPAR',thetname)
5715 open (ithep,file=thetname,status='old',action='read')
5716 ! print *,"Processor",myrank," opened file ITHEP"
5717 call getenv_loc('ROTPAR',rotname)
5718 open (irotam,file=rotname,status='old',action='read')
5720 call getenv_loc('ROTPAR_END',rotname_end)
5721 open (irotam_end,file=rotname_end,status='old',action='read')
5723 ! print *,"Processor",myrank," opened file IROTAM"
5724 call getenv_loc('TORPAR',torname)
5725 open (itorp,file=torname,status='old',action='read')
5726 ! print *,"Processor",myrank," opened file ITORP"
5727 call getenv_loc('TORDPAR',tordname)
5728 open (itordp,file=tordname,status='old',action='read')
5729 ! print *,"Processor",myrank," opened file ITORDP"
5730 call getenv_loc('SCCORPAR',sccorname)
5731 open (isccor,file=sccorname,status='old',action='read')
5732 ! print *,"Processor",myrank," opened file ISCCOR"
5733 call getenv_loc('FOURIER',fouriername)
5734 open (ifourier,file=fouriername,status='old',action='read')
5735 ! print *,"Processor",myrank," opened file IFOURIER"
5736 call getenv_loc('ELEPAR',elename)
5737 open (ielep,file=elename,status='old',action='read')
5738 ! print *,"Processor",myrank," opened file IELEP"
5739 call getenv_loc('SIDEPAR',sidename)
5740 open (isidep,file=sidename,status='old',action='read')
5742 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5743 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5744 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5745 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5746 call getenv_loc('TORPAR_NUCL',torname_nucl)
5747 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5748 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5749 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5750 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5751 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5753 call getenv_loc('LIPTRANPAR',liptranname)
5754 open (iliptranpar,file=liptranname,status='old',action='read')
5755 call getenv_loc('TUBEPAR',tubename)
5756 open (itube,file=tubename,status='old',action='read')
5757 call getenv_loc('IONPAR',ionname)
5758 open (iion,file=ionname,status='old',action='read')
5759 call getenv_loc('IONPAR_TRAN',iontranname)
5760 open (iiontran,file=iontranname,status='old',action='read')
5761 ! print *,"Processor",myrank," opened file ISIDEP"
5762 ! print *,"Processor",myrank," opened parameter files"
5764 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5765 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5766 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5767 ! Get parameter filenames and open the parameter files.
5768 call getenv_loc('BONDPAR',bondname)
5769 open (ibond,file=bondname,status='old')
5770 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5771 open (ibond_nucl,file=bondname_nucl,status='old')
5773 call getenv_loc('THETPAR',thetname)
5774 open (ithep,file=thetname,status='old')
5775 call getenv_loc('ROTPAR',rotname)
5776 open (irotam,file=rotname,status='old')
5778 call getenv_loc('ROTPAR_END',rotname_end)
5779 open (irotam_end,file=rotname_end,status='old')
5781 call getenv_loc('TORPAR',torname)
5782 open (itorp,file=torname,status='old')
5783 call getenv_loc('TORDPAR',tordname)
5784 open (itordp,file=tordname,status='old')
5785 call getenv_loc('SCCORPAR',sccorname)
5786 open (isccor,file=sccorname,status='old')
5787 call getenv_loc('FOURIER',fouriername)
5788 open (ifourier,file=fouriername,status='old')
5789 call getenv_loc('ELEPAR',elename)
5790 open (ielep,file=elename,status='old')
5791 call getenv_loc('SIDEPAR',sidename)
5792 open (isidep,file=sidename,status='old')
5794 open (ithep_nucl,file=thetname_nucl,status='old')
5795 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5796 open (irotam_nucl,file=rotname_nucl,status='old')
5797 call getenv_loc('TORPAR_NUCL',torname_nucl)
5798 open (itorp_nucl,file=torname_nucl,status='old')
5799 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5800 open (itordp_nucl,file=tordname_nucl,status='old')
5801 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5802 open (isidep_nucl,file=sidename_nucl,status='old')
5804 call getenv_loc('LIPTRANPAR',liptranname)
5805 open (iliptranpar,file=liptranname,status='old')
5806 call getenv_loc('TUBEPAR',tubename)
5807 open (itube,file=tubename,status='old')
5808 call getenv_loc('IONPAR',ionname)
5809 open (iion,file=ionname,status='old')
5810 call getenv_loc('IONPAR_NUCL',ionnuclname)
5811 open (iionnucl,file=ionnuclname,status='old')
5812 call getenv_loc('IONPAR_TRAN',iontranname)
5813 open (iiontran,file=iontranname,status='old')
5814 call getenv_loc('WATWAT',iwaterwatername)
5815 open (iwaterwater,file=iwaterwatername,status='old')
5816 call getenv_loc('WATPROT',iwaterscname)
5817 open (iwatersc,file=iwaterscname,status='old')
5820 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5822 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5823 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5824 ! Get parameter filenames and open the parameter files.
5825 call getenv_loc('BONDPAR',bondname)
5826 open (ibond,file=bondname,status='old',action='read')
5827 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5828 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5829 call getenv_loc('THETPAR',thetname)
5830 open (ithep,file=thetname,status='old',action='read')
5831 call getenv_loc('ROTPAR',rotname)
5832 open (irotam,file=rotname,status='old',action='read')
5834 call getenv_loc('ROTPAR_END',rotname_end)
5835 open (irotam_end,file=rotname_end,status='old',action='read')
5837 call getenv_loc('TORPAR',torname)
5838 open (itorp,file=torname,status='old',action='read')
5839 call getenv_loc('TORDPAR',tordname)
5840 open (itordp,file=tordname,status='old',action='read')
5841 call getenv_loc('SCCORPAR',sccorname)
5842 open (isccor,file=sccorname,status='old',action='read')
5844 call getenv_loc('THETPARPDB',thetname_pdb)
5845 print *,"thetname_pdb ",thetname_pdb
5846 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5847 print *,ithep_pdb," opened"
5849 call getenv_loc('FOURIER',fouriername)
5850 open (ifourier,file=fouriername,status='old',readonly)
5851 call getenv_loc('ELEPAR',elename)
5852 open (ielep,file=elename,status='old',readonly)
5853 call getenv_loc('SIDEPAR',sidename)
5854 open (isidep,file=sidename,status='old',readonly)
5856 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5857 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5858 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5859 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5860 call getenv_loc('TORPAR_NUCL',torname_nucl)
5861 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5862 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5863 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5864 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5865 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5866 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5867 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5868 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5869 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5870 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5871 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5872 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5873 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5876 call getenv_loc('LIPTRANPAR',liptranname)
5877 open (iliptranpar,file=liptranname,status='old',action='read')
5878 call getenv_loc('LIPBOND',lipbondname)
5879 open (ilipbond,file=lipbondname,status='old',action='read')
5880 call getenv_loc('LIPNONBOND',lipnonbondname)
5881 open (ilipnonbond,file=lipnonbondname,status='old',action='read')
5882 call getenv_loc('LIPPROT',lipprotname)
5883 open (imartprot,file=lipprotname,status='old',action='read')
5885 call getenv_loc('TUBEPAR',tubename)
5886 open (itube,file=tubename,status='old',action='read')
5887 call getenv_loc('IONPAR',ionname)
5888 open (iion,file=ionname,status='old',action='read')
5889 call getenv_loc('IONPAR_NUCL',ionnuclname)
5890 open (iionnucl,file=ionnuclname,status='old',action='read')
5891 call getenv_loc('IONPAR_TRAN',iontranname)
5892 open (iiontran,file=iontranname,status='old',action='read')
5893 call getenv_loc('WATWAT',iwaterwatername)
5894 open (iwaterwater,file=iwaterwatername,status='old',action='read')
5895 call getenv_loc('WATPROT',iwaterscname)
5896 open (iwatersc,file=iwaterscname,status='old',action='read')
5899 call getenv_loc('ROTPARPDB',rotname_pdb)
5900 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5903 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5904 #if defined(WINIFL) || defined(WINPGI)
5905 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5906 #elif (defined CRAY) || (defined AIX)
5907 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5909 open (iscpp_nucl,file=scpname_nucl,status='old')
5911 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5916 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5917 ! Use -DOLDSCP to use hard-coded constants instead.
5919 call getenv_loc('SCPPAR',scpname)
5920 #if defined(WINIFL) || defined(WINPGI)
5921 open (iscpp,file=scpname,status='old',readonly,shared)
5922 #elif (defined CRAY) || (defined AIX)
5923 open (iscpp,file=scpname,status='old',action='read')
5925 open (iscpp,file=scpname,status='old')
5927 open (iscpp,file=scpname,status='old',action='read')
5930 call getenv_loc('PATTERN',patname)
5931 #if defined(WINIFL) || defined(WINPGI)
5932 open (icbase,file=patname,status='old',readonly,shared)
5933 #elif (defined CRAY) || (defined AIX)
5934 open (icbase,file=patname,status='old',action='read')
5936 open (icbase,file=patname,status='old')
5938 open (icbase,file=patname,status='old',action='read')
5941 ! Open output file only for CG processes
5942 ! print *,"Processor",myrank," fg_rank",fg_rank
5943 if (fg_rank.eq.0) then
5945 if (nodes.eq.1) then
5948 npos = dlog10(dfloat(nodes-1))+1
5950 if (npos.lt.3) npos=3
5951 write (liczba,'(i1)') npos
5952 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5954 write (liczba,form) me
5955 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5956 liczba(:ilen(liczba))
5957 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5959 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5961 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5962 liczba(:ilen(liczba))//'.mol2'
5963 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5964 liczba(:ilen(liczba))//'.stat'
5966 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5967 //liczba(:ilen(liczba))//'.stat')
5968 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5971 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5972 liczba(:ilen(liczba))//'.const'
5977 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5978 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5979 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5980 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5981 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5983 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5985 rest2name=prefix(:ilen(prefix))//'.rst'
5987 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5990 #if defined(AIX) || defined(PGI)
5991 if (me.eq.king .or. .not. out1file) &
5992 open(iout,file=outname,status='unknown')
5994 if (fg_rank.gt.0) then
5995 write (liczba,'(i3.3)') myrank/nfgtasks
5996 write (ll,'(bz,i3.3)') fg_rank
5997 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
6002 open(igeom,file=intname,status='unknown',position='append')
6003 open(ipdb,file=pdbname,status='unknown')
6004 open(imol2,file=mol2name,status='unknown')
6005 open(istat,file=statname,status='unknown',position='append')
6007 !1out open(iout,file=outname,status='unknown')
6010 if (me.eq.king .or. .not.out1file) &
6011 open(iout,file=outname,status='unknown')
6013 if (fg_rank.gt.0) then
6014 write (liczba,'(i3.3)') myrank/nfgtasks
6015 write (ll,'(bz,i3.3)') fg_rank
6016 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
6021 open(igeom,file=intname,status='unknown',access='append')
6022 open(ipdb,file=pdbname,status='unknown')
6023 open(imol2,file=mol2name,status='unknown')
6024 open(istat,file=statname,status='unknown',access='append')
6026 !1out open(iout,file=outname,status='unknown')
6029 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
6030 csa_seed=prefix(:lenpre)//'.CSA.seed'
6031 csa_history=prefix(:lenpre)//'.CSA.history'
6032 csa_bank=prefix(:lenpre)//'.CSA.bank'
6033 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
6034 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
6035 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
6036 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
6037 csa_int=prefix(:lenpre)//'.int'
6038 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
6039 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
6040 csa_in=prefix(:lenpre)//'.CSA.in'
6041 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
6044 write (iout,'(80(1h-))')
6045 write (iout,'(30x,a)') "FILE ASSIGNMENT"
6046 write (iout,'(80(1h-))')
6047 write (iout,*) "Input file : ",&
6048 pref_orig(:ilen(pref_orig))//'.inp'
6049 write (iout,*) "Output file : ",&
6050 outname(:ilen(outname))
6052 write (iout,*) "Sidechain potential file : ",&
6053 sidename(:ilen(sidename))
6055 write (iout,*) "SCp potential file : ",&
6056 scpname(:ilen(scpname))
6058 write (iout,*) "Electrostatic potential file : ",&
6059 elename(:ilen(elename))
6060 write (iout,*) "Cumulant coefficient file : ",&
6061 fouriername(:ilen(fouriername))
6062 write (iout,*) "Torsional parameter file : ",&
6063 torname(:ilen(torname))
6064 write (iout,*) "Double torsional parameter file : ",&
6065 tordname(:ilen(tordname))
6066 write (iout,*) "SCCOR parameter file : ",&
6067 sccorname(:ilen(sccorname))
6068 write (iout,*) "Bond & inertia constant file : ",&
6069 bondname(:ilen(bondname))
6070 write (iout,*) "Bending parameter file : ",&
6071 thetname(:ilen(thetname))
6072 write (iout,*) "Rotamer parameter file : ",&
6073 rotname(:ilen(rotname))
6076 write (iout,*) "Thetpdb parameter file : ",&
6077 thetname_pdb(:ilen(thetname_pdb))
6080 write (iout,*) "Threading database : ",&
6081 patname(:ilen(patname))
6083 write (iout,*)" DIRTMP : ",&
6085 write (iout,'(80(1h-))')
6088 end subroutine openunits
6089 !-----------------------------------------------------------------------------
6092 use geometry_data, only: nres,dc
6094 ! implicit real*8 (a-h,o-z)
6095 ! include 'DIMENSIONS'
6096 ! include 'COMMON.CHAIN'
6097 ! include 'COMMON.IOUNITS'
6098 ! include 'COMMON.MD'
6101 ! real(kind=8) :: var,ene
6103 open(irest2,file=rest2name,status='unknown')
6104 read(irest2,*) totT,EK,potE,totE,t_bath
6107 ! AL 4/17/17: Now reading d_t(0,:) too
6109 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
6112 ! AL 4/17/17: Now reading d_c(0,:) too
6114 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
6117 read (irest2,*) iset
6121 end subroutine readrst
6122 !-----------------------------------------------------------------------------
6123 subroutine read_fragments
6127 use control_data, only:out1file
6130 ! implicit real*8 (a-h,o-z)
6131 ! include 'DIMENSIONS'
6135 ! include 'COMMON.SETUP'
6136 ! include 'COMMON.CHAIN'
6137 ! include 'COMMON.IOUNITS'
6138 ! include 'COMMON.MD'
6139 ! include 'COMMON.CONTROL'
6142 ! real(kind=8) :: var,ene
6144 read(inp,*) nset,nfrag,npair,nfrag_back
6146 !el from module energy
6147 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
6148 if(.not.allocated(wfrag_back)) then
6149 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6150 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
6152 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
6153 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
6155 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
6156 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
6159 if(me.eq.king.or..not.out1file) &
6160 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
6161 " nfrag_back",nfrag_back
6163 read(inp,*) mset(iset)
6165 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
6167 if(me.eq.king.or..not.out1file) &
6168 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
6169 ifrag(2,i,iset), qinfrag(i,iset)
6172 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
6174 if(me.eq.king.or..not.out1file) &
6175 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
6176 ipair(2,i,iset), qinpair(i,iset)
6179 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6180 wfrag_back(3,i,iset),&
6181 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6182 if(me.eq.king.or..not.out1file) &
6183 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
6184 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
6188 end subroutine read_fragments
6189 !-----------------------------------------------------------------------------
6191 !-----------------------------------------------------------------------------
6195 ! implicit real*8 (a-h,o-z)
6196 ! include 'DIMENSIONS'
6197 ! include 'COMMON.CSA'
6198 ! include 'COMMON.BANK'
6199 ! include 'COMMON.IOUNITS'
6201 ! integer :: ntf,ik,iw_pdb
6202 ! real(kind=8) :: var,ene
6204 open(icsa_in,file=csa_in,status="old",err=100)
6205 read(icsa_in,*) nconf
6206 read(icsa_in,*) jstart,jend
6207 read(icsa_in,*) nstmax
6208 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6209 read(icsa_in,*) nran0,nran1,irr
6210 read(icsa_in,*) nseed
6211 read(icsa_in,*) ntotal,cut1,cut2
6212 read(icsa_in,*) estop
6213 read(icsa_in,*) icmax,irestart
6214 read(icsa_in,*) ntbankm,dele,difcut
6215 read(icsa_in,*) iref,rmscut,pnccut
6216 read(icsa_in,*) ndiff
6223 end subroutine csa_read
6224 !-----------------------------------------------------------------------------
6225 subroutine initial_write
6228 ! implicit real*8 (a-h,o-z)
6229 ! include 'DIMENSIONS'
6230 ! include 'COMMON.CSA'
6231 ! include 'COMMON.BANK'
6232 ! include 'COMMON.IOUNITS'
6234 ! integer :: ntf,ik,iw_pdb
6235 ! real(kind=8) :: var,ene
6237 open(icsa_seed,file=csa_seed,status="unknown")
6238 write(icsa_seed,*) "seed"
6240 #if defined(AIX) || defined(PGI)
6241 open(icsa_history,file=csa_history,status="unknown",&
6244 open(icsa_history,file=csa_history,status="unknown",&
6247 write(icsa_history,*) nconf
6248 write(icsa_history,*) jstart,jend
6249 write(icsa_history,*) nstmax
6250 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6251 write(icsa_history,*) nran0,nran1,irr
6252 write(icsa_history,*) nseed
6253 write(icsa_history,*) ntotal,cut1,cut2
6254 write(icsa_history,*) estop
6255 write(icsa_history,*) icmax,irestart
6256 write(icsa_history,*) ntbankm,dele,difcut
6257 write(icsa_history,*) iref,rmscut,pnccut
6258 write(icsa_history,*) ndiff
6260 write(icsa_history,*)
6263 open(icsa_bank1,file=csa_bank1,status="unknown")
6264 write(icsa_bank1,*) 0
6268 end subroutine initial_write
6269 !-----------------------------------------------------------------------------
6270 subroutine restart_write
6273 ! implicit real*8 (a-h,o-z)
6274 ! include 'DIMENSIONS'
6275 ! include 'COMMON.IOUNITS'
6276 ! include 'COMMON.CSA'
6277 ! include 'COMMON.BANK'
6279 ! integer :: ntf,ik,iw_pdb
6280 ! real(kind=8) :: var,ene
6282 #if defined(AIX) || defined(PGI)
6283 open(icsa_history,file=csa_history,position="append")
6285 open(icsa_history,file=csa_history,access="append")
6287 write(icsa_history,*)
6288 write(icsa_history,*) "This is restart"
6289 write(icsa_history,*)
6290 write(icsa_history,*) nconf
6291 write(icsa_history,*) jstart,jend
6292 write(icsa_history,*) nstmax
6293 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
6294 write(icsa_history,*) nran0,nran1,irr
6295 write(icsa_history,*) nseed
6296 write(icsa_history,*) ntotal,cut1,cut2
6297 write(icsa_history,*) estop
6298 write(icsa_history,*) icmax,irestart
6299 write(icsa_history,*) ntbankm,dele,difcut
6300 write(icsa_history,*) iref,rmscut,pnccut
6301 write(icsa_history,*) ndiff
6302 write(icsa_history,*)
6303 write(icsa_history,*) "irestart is: ", irestart
6305 write(icsa_history,*)
6309 end subroutine restart_write
6310 !-----------------------------------------------------------------------------
6312 !-----------------------------------------------------------------------------
6313 subroutine write_pdb(npdb,titelloc,ee)
6315 ! implicit real*8 (a-h,o-z)
6316 ! include 'DIMENSIONS'
6317 ! include 'COMMON.IOUNITS'
6318 character(len=50) :: titelloc1
6319 character*(*) :: titelloc
6320 character(len=3) :: zahl
6321 character(len=5) :: liczba5
6323 integer :: npdb !,ilen
6327 ! real(kind=8) :: var,ene
6331 if (npdb.lt.1000) then
6332 call numstr(npdb,zahl)
6333 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
6335 if (npdb.lt.10000) then
6336 write(liczba5,'(i1,i4)') 0,npdb
6338 write(liczba5,'(i5)') npdb
6340 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
6342 call pdbout(ee,titelloc1,ipdb)
6345 end subroutine write_pdb
6346 !-----------------------------------------------------------------------------
6348 !-----------------------------------------------------------------------------
6349 subroutine write_thread_summary
6350 ! Thread the sequence through a database of known structures
6351 use control_data, only: refstr
6353 use energy_data, only: n_ene_comp
6355 ! implicit real*8 (a-h,o-z)
6356 ! include 'DIMENSIONS'
6358 use MPI_data !include 'COMMON.INFO'
6361 ! include 'COMMON.CONTROL'
6362 ! include 'COMMON.CHAIN'
6363 ! include 'COMMON.DBASE'
6364 ! include 'COMMON.INTERACT'
6365 ! include 'COMMON.VAR'
6366 ! include 'COMMON.THREAD'
6367 ! include 'COMMON.FFIELD'
6368 ! include 'COMMON.SBRIDGE'
6369 ! include 'COMMON.HEADER'
6370 ! include 'COMMON.NAMES'
6371 ! include 'COMMON.IOUNITS'
6372 ! include 'COMMON.TIME1'
6374 integer,dimension(maxthread) :: ip
6375 real(kind=8),dimension(0:n_ene) :: energia
6377 integer :: i,j,ii,jj,ipj,ik,kk,ist
6378 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
6380 write (iout,'(30x,a/)') &
6381 ' *********** Summary threading statistics ************'
6382 write (iout,'(a)') 'Initial energies:'
6383 write (iout,'(a4,2x,a12,14a14,3a8)') &
6384 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
6385 'RMSnat','NatCONT','NNCONT','RMS'
6386 ! Energy sort patterns
6391 enet=ener(n_ene-1,ip(i))
6394 if (ener(n_ene-1,ip(j)).lt.enet) then
6396 enet=ener(n_ene-1,ip(j))
6408 ist=nres_base(2,ii)+ipatt(2,i)
6410 energia(i)=ener0(kk,i)
6412 etot=ener0(n_ene_comp+1,i)
6413 rmsnat=ener0(n_ene_comp+2,i)
6414 rms=ener0(n_ene_comp+3,i)
6415 frac=ener0(n_ene_comp+4,i)
6416 frac_nn=ener0(n_ene_comp+5,i)
6419 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6420 i,str_nam(ii),ist+1,&
6421 (energia(print_order(kk)),kk=1,nprint_ene),&
6422 etot,rmsnat,frac,frac_nn,rms
6424 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
6425 i,str_nam(ii),ist+1,&
6426 (energia(print_order(kk)),kk=1,nprint_ene),etot
6429 write (iout,'(//a)') 'Final energies:'
6430 write (iout,'(a4,2x,a12,17a14,3a8)') &
6431 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
6432 'RMSnat','NatCONT','NNCONT','RMS'
6436 ist=nres_base(2,ii)+ipatt(2,i)
6438 energia(kk)=ener(kk,ik)
6440 etot=ener(n_ene_comp+1,i)
6441 rmsnat=ener(n_ene_comp+2,i)
6442 rms=ener(n_ene_comp+3,i)
6443 frac=ener(n_ene_comp+4,i)
6444 frac_nn=ener(n_ene_comp+5,i)
6445 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6446 i,str_nam(ii),ist+1,&
6447 (energia(print_order(kk)),kk=1,nprint_ene),&
6448 etot,rmsnat,frac,frac_nn,rms
6450 write (iout,'(/a/)') 'IEXAM array:'
6451 write (iout,'(i5)') nexcl
6453 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
6455 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
6456 'Max. time for threading step ',max_time_for_thread,&
6457 'Average time for threading step: ',ave_time_for_thread
6459 end subroutine write_thread_summary
6460 !-----------------------------------------------------------------------------
6461 subroutine write_stat_thread(ithread,ipattern,ist)
6463 use energy_data, only: n_ene_comp
6465 ! implicit real*8 (a-h,o-z)
6466 ! include "DIMENSIONS"
6467 ! include "COMMON.CONTROL"
6468 ! include "COMMON.IOUNITS"
6469 ! include "COMMON.THREAD"
6470 ! include "COMMON.FFIELD"
6471 ! include "COMMON.DBASE"
6472 ! include "COMMON.NAMES"
6473 real(kind=8),dimension(0:n_ene) :: energia
6475 integer :: ithread,ipattern,ist,i
6476 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
6478 #if defined(AIX) || defined(PGI)
6479 open(istat,file=statname,position='append')
6481 open(istat,file=statname,access='append')
6484 energia(i)=ener(i,ithread)
6486 etot=ener(n_ene_comp+1,ithread)
6487 rmsnat=ener(n_ene_comp+2,ithread)
6488 rms=ener(n_ene_comp+3,ithread)
6489 frac=ener(n_ene_comp+4,ithread)
6490 frac_nn=ener(n_ene_comp+5,ithread)
6491 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
6492 ithread,str_nam(ipattern),ist+1,&
6493 (energia(print_order(i)),i=1,nprint_ene),&
6494 etot,rmsnat,frac,frac_nn,rms
6497 end subroutine write_stat_thread
6498 !-----------------------------------------------------------------------------
6500 subroutine readpdb_template(k)
6501 ! Read the PDB file for read_constr_homology with read2sigma
6502 ! and convert the peptide geometry into virtual-chain geometry.
6504 ! include 'DIMENSIONS'
6505 ! include 'COMMON.LOCAL'
6506 ! include 'COMMON.VAR'
6507 ! include 'COMMON.CHAIN'
6508 ! include 'COMMON.INTERACT'
6509 ! include 'COMMON.IOUNITS'
6510 ! include 'COMMON.GEO'
6511 ! include 'COMMON.NAMES'
6512 ! include 'COMMON.CONTROL'
6513 ! include 'COMMON.FRAG'
6514 ! include 'COMMON.SETUP'
6515 use compare_data, only:nhfrag,nbfrag
6516 integer :: i,j,k,ibeg,ishift1,ires,iii,ires_old,ishift,ity, &
6518 logical lprn /.false./,fail
6519 real(kind=8), dimension (3):: e1,e2,e3
6520 real(kind=8) :: dcj,efree_temp
6524 real(kind=8), dimension (3,40) :: sccor
6526 integer, dimension (:), allocatable :: iterter
6527 if(.not.allocated(iterter))allocate(iterter(nres))
6534 write (2,*) "UNRES_PDB",unres_pdb
6542 read (ipdbin,'(a80)',end=10) card
6543 if (card(:3).eq.'END') then
6545 else if (card(:3).eq.'TER') then
6548 itype(ires_old-1,1)=ntyp1
6549 iterter(ires_old-1)=1
6550 #if defined(WHAM_RUN) || defined(CLUSTER)
6551 if (ires_old.lt.nres) then
6553 itype(ires_old,1)=ntyp1
6555 #if defined(WHAM_RUN) || defined(CLUSTER)
6559 ! write (iout,*) "Chain ended",ires,ishift,ires_old
6562 dc(j,ires)=sccor(j,iii)
6565 call sccenter(ires,iii,sccor)
6568 ! Fish out the ATOM cards.
6569 if (index(card(1:4),'ATOM').gt.0) then
6570 read (card(12:16),*) atom
6571 ! write (iout,*) "! ",atom," !",ires
6572 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
6573 read (card(23:26),*) ires
6574 read (card(18:20),'(a3)') res
6575 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
6576 ! & " ires_old",ires_old
6577 ! write (iout,*) "ishift",ishift," ishift1",ishift1
6578 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
6579 if (ires-ishift+ishift1.ne.ires_old) then
6580 ! Calculate the CM of the preceding residue.
6584 dc(j,ires_old)=sccor(j,iii)
6587 call sccenter(ires_old,iii,sccor)
6591 ! Start new residue.
6592 if (res.eq.'Cl-' .or. res.eq.'Na+') then
6595 else if (ibeg.eq.1) then
6596 ! write (iout,*) "BEG ires",ires
6598 if (res.ne.'GLY' .and. res.ne. 'ACE') then
6602 ires=ires-ishift+ishift1
6604 ! write (iout,*) "ishift",ishift," ires",ires,
6605 ! & " ires_old",ires_old
6606 ! write (iout,*) "ires",ires," ibeg",ibeg," ishift",ishift
6608 else if (ibeg.eq.2) then
6610 ishift=-ires_old+ires-1
6612 ! write (iout,*) "New chain started",ires,ishift
6615 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
6616 ires=ires-ishift+ishift1
6619 if (res.eq.'ACE' .or. res.eq.'NHE') then
6622 itype(ires,1)=rescode(ires,res,0,1)
6625 ires=ires-ishift+ishift1
6627 ! write (iout,*) "ires_old",ires_old," ires",ires
6628 ! if (card(27:27).eq."A" .or. card(27:27).eq."B") then
6631 ! write (2,*) "ires",ires," res ",res," ity",ity
6632 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
6633 res.eq.'NHE'.and.atom(:2).eq.'HN') then
6634 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
6635 ! write (iout,*) "backbone ",atom ,ires,res, (c(j,ires),j=1,3)
6637 write (iout,'(2i3,2x,a,3f8.3)') &
6638 ires,itype(ires,1),res,(c(j,ires),j=1,3)
6642 sccor(j,iii)=c(j,ires)
6644 if (ishift.ne.0) then
6645 ires_ca=ires+ishift-ishift1
6649 ! write (*,*) card(23:27),ires,itype(ires)
6650 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and.&
6651 atom.ne.'N' .and. atom.ne.'C' .and.&
6652 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and.&
6653 atom.ne.'OXT' .and. atom(:2).ne.'3H') then
6654 ! write (iout,*) "sidechain ",atom
6656 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
6660 10 if(me.eq.king.or..not.out1file) &
6661 write (iout,'(a,i5)') ' Nres: ',ires
6662 ! Calculate dummy residue coordinates inside the "chain" of a multichain
6665 write(2,*) "tutaj",ires,nres
6667 ! write (iout,*) i,itype(i),itype(i+1)
6668 if (itype(i,1).eq.ntyp1.and.iterter(i).eq.1) then
6669 if (itype(i+1,1).eq.ntyp1.and.iterter(i+1).eq.1 ) then
6670 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
6671 ! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
6672 ! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
6674 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6675 call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
6682 c(j,i)=c(j,i-1)-1.9d0*e2(j)
6686 dcj=(c(j,i-2)-c(j,i-3))/2.0
6687 if (dcj.eq.0) dcj=1.23591524223
6692 else !itype(i+1).eq.ntyp1
6694 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6695 call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
6702 c(j,i)=c(j,i+1)-1.9d0*e2(j)
6706 dcj=(c(j,i+3)-c(j,i+2))/2.0
6707 if (dcj.eq.0) dcj=1.23591524223
6712 endif !itype(i+1).eq.ntyp1
6713 endif !itype.eq.ntyp1
6715 ! Calculate the CM of the last side chain.
6718 dc(j,ires)=sccor(j,iii)
6721 call sccenter(ires,iii,sccor)
6725 if (itype(nres,1).ne.10) then
6729 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
6730 call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
6737 c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
6741 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
6742 if (dcj.eq.0) dcj=1.23591524223
6743 c(j,nres)=c(j,nres-1)+dcj
6744 c(j,2*nres)=c(j,nres)
6755 c(j,2*nres)=c(j,nres)
6757 if (itype(1,1).eq.ntyp1) then
6761 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
6762 call refsys(2,3,4,e1,e2,e3,fail)
6769 c(j,1)=c(j,2)-1.9d0*e2(j)
6773 dcj=(c(j,4)-c(j,3))/2.0
6779 ! Copy the coordinates to reference coordinates
6785 ! Calculate internal coordinates.
6786 if (out_template_coord) then
6787 write (iout,'(/a)') &
6788 "Cartesian coordinates of the reference structure"
6789 write (iout,'(a,3(3x,a5),5x,3(3x,a5))') &
6790 "Residue ","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
6792 write (iout,'(a3,1x,i4,3f8.3,5x,3f8.3)')&
6793 restyp(itype(ires,1),1),ires,(c(j,ires),j=1,3),&
6794 (c(j,ires+nres),j=1,3)
6797 ! Calculate internal coordinates.
6798 call int_from_cart(.true.,out_template_coord)
6799 call sc_loc_geom(.false.)
6801 thetaref(i)=theta(i)
6806 dc(j,i)=c(j,i+1)-c(j,i)
6807 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
6812 dc(j,i+nres)=c(j,i+nres)-c(j,i)
6813 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
6815 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),
6816 ! & vbld_inv(i+nres)
6821 cref(j,i+nres,1)=c(j,i+nres)
6831 end subroutine readpdb_template
6832 !-----------------------------------------------------------------------------
6834 !-----------------------------------------------------------------------------
6835 end module io_config